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FOREWORD 


This is a combined report for the two project* at the Massachusetts Institute of Technology which are 
uponuored by the Office of Naval Rraearch under Contract N5orl60. 


PART 1 

Machine Methods of Computation and Numerical Analysis 
1. GENERAL COMMENTS 


Project on Machine Methods of Computation and Numerical Analysis 

Thia Project is an outgrowth of the activities of the Institute Committee on Machine Methods of Compute* 
tlon, established in November 1950. The purpose of the Project is (1) to integrate the efforts of all the 
departments and groups at MIT who are working with modern computing machines and their applications, and 
(2) to train men in the use of these machines for computation and numerical analysis. 

People from several departments of the Institute are taking part in the project. In the Appendix will be 
found a list of the personnel active in this program. 


Project Whirlwind 

This Project makes use of the facilities of the Digital Computer Laboratory. The principal objective of 
the Project is the application of an electronic digital computer of large capacity and very high speed (Whirl* 
wind 1) to problems in mathematics, science, engineering, stimulation, and control. 


The Whirlwind 1 Computer 

Whirlwind I is of the high-speed electronic digital type, in which quantities are represented as discrete 
numbers, and complex problems are solved by the repeated use of fundamental arithmetic and logical (i.e., 
control or selection) operations. Computations are executed by fractional-microsecond pulses in electronic 
circuits, of which the principal ones are U) the flip-flop, a circuit containing two vacuum tubes so connected 
that one tube or the other is conducting, but not both; (2) the gate or coincidence circuit; (3) the magnetic-core 
memory, in which binary digits are stored as one of two directions of magnetic flux within ferro-magnetic 
cores. 

Whirlwind i uses numbers of 16 binary digits (equivalent to about 5 decimal digits). This length was 
selected to limit the machine to a practical size, but it permits the computation of many simulation problems. 
Calculations requiring greater number length are handled by the use of multiple-length numbers. Rapid-access 
magnetic-core memory has a capacity of 32,768 binary digits. Present speed of the computer is 40,000 Bingle- 
address operalibns per second, equivalent to about 20.000 multiplications per second. This speed is higher than 
general scientific computation demands at the present state of the art. but is needed for control and simulation 
studies. 


During the fail and winter the contacts between the project staff and others workuig in the field of 
numerical analysis and machine computation have increased in frequency and broadened in scope. The Discus¬ 
sion Group on Numerical Analysis to be mentioned in section 3, has stimulated the interest and participation of 
a number of members of the Mathematics Department. Other research projects finding a need for machine 
calculation, huve come to depend on thia project for help in finding how best to use computing equipment. The 
Solid-State and Molecular Theory Group under Professor J. C. Slater, has been using Whirlwind and other 
computing equipment to an increasing extont. Corbato, a member of the Machine Project, has been spending 
much of hia time with Slater 1 a Group, helping various members code their problems and working out appro¬ 
priate computing methods. Uretsky, Sartori and Rotenberg have been working with members of the Theoretical 
Group of the Research Laboratory for Nuclear Science and Engineering, applying machine techniques to numerous 
problems in this field. 

There has also been close contact with the staff of the Operations Research Project, particularly in the 
study of problems of Linear Programming and techniques for their solution. Joint seminars have been held. 

Some of the Whirlwind staff have been examining the simplex method and other techniques for obtaining solutions 
to linear programming problems. John D. C. Little, who transferred from this project to the Operations Research 
project last year, completed the solution to a dynamic programming problem this fall, using Whirlwind to obtain 
numerical results. Other groups, both in engineering and in pure science, have received advice and help in 
connection with their computing or data handling problems. 

The research projects reported in the following pages, represent a wide variety of subjects; their 
common aspects lie in the adaptation of the problem for machine computation. Each application has added to 
our experience in the techniques of such application, many of them have resulted in useful sub-routines, which 
arc now In the Whirlwind library, available to other workers. 



4 


5 








APPROVED FOR PUBLIC RELEASE. CASE 06-1104 


2. GRADUATE SCHOOL RESEARCH 


2.1 Index to Reports 

Title Page 

A Study of the Basic Problem of Numcrtal Analysis Expressed ui the Language of 

Computing Machines 7 

Evaluation of the Explicit Difference Formula for a Parabolic Differential Equation H 

Calculation of number of Structures of Relations on Finite Sets IU 

The Stability of Thin, Shallow F.lastic Shells 12 

Machine Solution of the Diffusion Equation 13 

Multi-center Integrals | u 

Eigen values for a Spheroidal Square Well 10 

Variational Determination of Atomic Wave Functions 21 

Neutron-Deuteron Scattering 23 

Nuclear Constitution 29 

Coulomb Wave Functions 30 

Computation of Chemical Engineering Problems of Multistage Distillation, Absorption, 31 

and Extraction 

Response of a Five-Story Frame Building to Dynamic Loading 31 

Analysis of a T>o Story Steel Frame Building Under Dynamic Loading 31 

Reactor Runaway Prevention 33 

Newton Analysis of Earth Resistivity Measurements 33 


GRADUATE SCHOOL RESEARCH 


2.2 Progress Reports 

A STUDY OF THE BASIC PROBLEM OF NUMERICAL ANALYSIS EXPRESSED IN THE LANGUAGE OF 

COMPUTING MACHINES 

The following report introduces a new project that has developed partly out of the discussion group on 
numerical analysis referred to under Group Activities in the last quarterly report. The discussion group is 
still continuing (as recorded in section 3 of the present quarterly) and will be used as a place to work over and 
extend the ideas that are presented here. 

The question of what aspects of numerical analysis are most important for improving machine methods 
is not much different from the question of what are the essential aspects of numerical analysis. The reason for 
this is that numerical analysis is to be regarded as a study of computational methods, or methods for solving 
problems that ultimately yield numbers, and as such is directed toward computing machines. Of course, the 
mathematn lan need not distinguish between a calculutional procedure that is carried out by hand and one that 
is carried out mechanically or electronically; both are procedures for a machine. Therefore, it seems only 
realistti thut numerical analysis should be studied in conjunction with the logic of computing machines and that 
the problems of numerical analysis should be couched in the language of the computing machines for which they 
are destined. 

This observation should be balanced with the realization that the language of mathematics as opposed 
to the language of machines, is not highly adapted to computational procedures except in the most simple 
cases of arithmetic. The Language of mathematics in which most numerical analysis is done accomodates 
concepts of limits, for example, which are present in computations only in an approximate sense. On the other 
hand, the concept of time required to carry out a mathematical operation, or the concept of deciuion (so essen¬ 
tial to iteration procedures) are only indirectly accomodated by the usual tools of mathematics and at best 
carry a very different significance in the fields of abstract and applied mathematics. 

Thus, while numerical analysis has an affinity for computing machines, the language in which it is studied 
is in some ways foreign to them. There are two unfortunate results: first, the basic problems of numerical 
analysis arc at times obscured or complicated, and, second, all results must be translated into the language of 
machines. In the present paper there is an attempt to rectify this situation by defining in the language of 
machines the most basic problem of numerical analysis. We will later attempt to outline some techniques for 
solving the problem. Such an outline can only be tentative and any ideas in it can at best be exemplary. It should 
also be emphasized that the "basic problem of numerical analysis" is to some extent born of a subjective 
judgement. 

In the language of computing machines, a numerical solution to a mathematical problem is a collection 
of parameter values and a sequence of Instructions. The parameters and instructions determine an output in 
the form of a set of numbers. The numerical solution to a particular problem solved on a particular machine 
lias, of course, the same ingredients, but the instructions must be specialized to the machine and the output 
must be close in some sense to a desired set of numbers. (It is unrealistic to speak of exactness of an output 
to a desired set of numbers, because only in rare cases can a numerical solution be an exact solution to a 
problem.) Having defined a numerical solution, it is clear that there are many such for each mathematical 
problem, because there arc many sequences of instructions with outputs that are within a certain degree of 
closeness to a desired set of numbers. The flexibility that permits a multiplicity of numerical solutions comes, 
first, from the fact that many sequences of instructions give the same output and, second, from the fact that 
the same output is not even required in order for two numerical solutions to satisfy a given problem. Thus if 
one attaches a measure of time to each sequence (the length of time for the machine to perform the sequence) 
and a measure of accuracy (the closeness of the output to a desired set of values) there will be two quantitative 
ways to compare the many numerical solutions to a given problem. 

It is proposed here that numerical analysis is not chiefly concerned with finding a possible numerical 
solution among a great multiplicity. The basic problem of numeric al analysis is to find the fastest and most 
accurate out of a set of possible numerical solutions. This statement will now be made more precise. 

It is first necessary to define mathematically what in meant hv s machine. M. Let P - 
be a sequence chosen from the space, , of sequences of real valued parameters. P represents the input 

data or other parameters initially held in the memory of the machine. Let I «(Ij, l^.I ) he a sequence 

chosen from the space,, of sequences of Instructions. The instructions are assumed to pertain to the 
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mat him* M. Let O be a mapping of the pairs (P, I) onto the space of sequences (O^. O . . O^) of real numbers, 

O (P. I) is the output determined by the parameter set P and the instructions L Let T be a mapping of the pairs 
(P, 1) onto the positive real axis. T (P, 1) is the time required for the machine M to perform the sequence I, 
assuming that the parameters in the machine's memory are P. Now for mathematical purposes, the machine M 
is completely specified by the cross-product space and the two mappings defined on this space. O and T. 

Thus, we will write: 

M = ( 9* $• , O, T>. 


In the form of 


( 2 ) 




where r f is the ratio for the formula which may be different from the real time ratio r. The exact 

solution for equation 0) is in the form of 


The points (P, 1) in the space V>$ may be described as possible numerical solutions to a mathemati¬ 
cal problem. Of course, one solution may be better than another for a particular problem, and we ultimately 
want to say in what sense a possible numerical solution is optimal for a particular mathematical problem. 
Therefore, wc must now define an entity which characterises each mathematical problem for the machine M. 

This is determined by two components, a possible solution (P . 1 ) whose output 0(P , 1 ) is a desirable set of 
numbers for the mathematical problem, and a distance function, £, defined on the space of all possible pairs 
of outputs [ 0(P, II, O <P'. 1*) J , The entity, itself, is a real valued function. D, defined on the cross-product 
space 9*$ : 

DIP, D • $[ OIP, 0, OIP.ll] . 

o o 

It is necessary to assume that (P^, I ) is known in order to state ttie baste problem of numerical analysis 
in the language of machines. This assumption Is Justified, because numerical analysis is chiefly concerned with 
something deepet than finding a set of numbers that approximately satisfy the restrictions imposed by a mathe¬ 
matical problem. After all, it la not difficult, if there ore no practical limitations, to specify a sequence of 
instructions which yield numbers closely approximating an integral of a function, the derivative of a function, the 
solution of an integral equation, or the solution of a differential equation. Thus. (P . 1 ) may be viewed as the 
solution to a mathematical problem translated into machine language without any regard for practicality, auch 
as Speed of convergence, number of iterations Involved, etc.. The distance function may be viewed as a criterion 
for Judging how a solution that may be more practical than (P , 1 ) compares in accuracy to (P I ). 

O O * o’ o 

After these considerations, the basic problem of numerical analysis appears to have the following two 

forms: 


Form I: Given the machine < § , O, T), the mathematical problem D, and the positive number D , 

to find (P*. I') such that TIP'. I') • min T(P, 1 ). where S_^ • J (p # D: D (P. I) < D > ° 

D o ‘ ° 1 

FormJDh Given the machine ( . O, T), the mathematical problem D. and the positive number T 

to find (P . 1 ) such that D(P’\ I") * mm D (P. I), where S . ( (P, 0: T(P, I) < T f . ° 

a T o 4 o * 


T. - w * t / _ » 

(3) t « v; ^ ♦ £ c k e -* ***> 

• m* 

and that for equation (2) is 


where X K - l - th( l - 61 
k - r/ar ' 

/ * 

and O and are characteristic values depending on the boundary conditions. 

A method to optimise the explicit approximation for coarse grid network Ithst is, M * J/aals small) is explained 
in the Quarterly Progress Report No. 13 (September 15. 1954). Initial conditions are , - 0 • Ths 

boundary conditions are } y, - o 


and (1) y, t ' 1 

<“> 

un> vi e - » 

The optimized parameters ry r. and VZ, or Y, 0 
Progress Report No. 14 (January 1, 1955). 


where Q * 1, and 3. 

the corner values for cases 1 and U are given in Quarterly 


For case III. the second order difference form for the boundary condition is 

< 5 > Vl.a - *(>> ^ 1 • 

The PJs in (4) have all the M*1 values instead of M. The characteristic values are the roots of the trams 
eendcntal equation 


If the mmimums do not exist or are not unique, of course the above rorms must be slightly modified. 

It should be pointed out that in any practical case, the outputs CXP, l), OIP , 1 1 will not be known to the 
numerical analyst when he attempts his problem of optimizing. Nevertheless, it while shown later how DIP I) 
can be determined once j, and <P I ) are specified and how the minimisation can be carried out at least ’ 

; h 7,:“ y n U “ l ““ Um * T(P - l> •* ««»>«>' ordering relations of 

the type TIP I) < TIP I ). will in general be known. In fact, though the above forma of the basic problem 

aSTn“n^rto^X“ ' “ U knowled * e l ‘ bou, ° and T ‘ h <“ ^ntifiea the problem with machines 


Bayard Rankin 

EVALUATION OF THE EXPLICIT DIFFERENCE FORMULA FOR A PARABOLIC DIFFERENTIAL EQUATION 
The present report is concerned with the explicit difference approxim.tlon to the diffusion equation 

(1) iLj? * 2-* , O t a i I a i -t 

Sx‘ Jt > 


(6) /“\(?^ f*JK S "P- m * • 

The M+l th value is fS M<J . n> z 1 y which is complex. Similarly, the rfn'sln (3) are the roots of the trans- 
eendcntal equation 

(7) (2 « k r,„ » 1 

Since at n' a from (7) and 1 s from (6) have no simple relationship, direct comparison of (3) and (4) shows 
that not only the excitation modes. C k and , and the decaying modes. C ”‘ r and ,1* are different but 
also the space modes are different. The optimization scheme is to make 

(1) D ( » C ( by choosing 

(21 Tt k a c ‘ *" ky/ng for n - 1 and 2 by choosing r f and r 

and (3) • «■<, by choosing Q f different from Q in equations (6) and (7). 

The merits of this optimization method are evaluated by comparing results of the same condtUer •.fltb 
r . i [ l]and with the evret Uto^Ucal solution. Step by step calculations are made by the Whirlwind Computer. 

S 








APPROVED FOR PUBLIC RELEASE. CASE 06-1104. 


GRADUATE SCHOOL RESEARCH 


GRADUATE SCHOOL RESEARCH 


Results of the rases I and U fco far show that the excitation modes are more important when r is in 
the neighborhood of ^ . Both the present method and the r « g with a corner value adjustment yield 
simiUirly good results and are both superior to cases where r is < housen to be other values. For case 111 with 
Q * - . the present method gives many times better results. 

Results are still being completed and studied. Investigation also includes application of the step responses 
to varying boundary conditions, such us t ■ 5i« l t , by the method of superposition. Details of the 
complete results will be available in the filial report in the corning Quarterly. 

Stubility . rite non can be shown by equation (4). the explicit formula will decay toward the steady state 
value when J j i 1. In other words 

(8) Vi I/O - . 

The critical value lies on the highest mode. ThuH, for case 1 where y is prescribed. 

2W> TP ant * r - 1 


The possible numberings of the set an- simply the n! permutations forming the symmetric group of 
order nl Consider the matrix transformation which occurs when surh ;* permutation is effected on members of 
the set without altering the relational structure. Let the permutation be represented by JT. renumbering obje. t 
l with the label « • TM. Then if the relational matrix is A. this permutation introduces a transformation, 
of the matrix A into a new matrix A' . Symbolically, 

r r A . a- 

“ij * rr, 

Thus uny basic matrix pattern muy be transformed into nl ostensibly different patterns by application of 
e.n h of the n! induced transformations 7«i . Of course some of these n! patterns may be repetitions of each 
other, for example the matrix all of whose elements are one transforms into itself under every . If we list 
all the 2 n ^ matrices and perform the n! transformations on each, we find the list partitioned into exclusive and 
exhaustive subgroups or "orbits". Every member matrix of one orbit may be obtained from every other member 
by applying a suitable 1^ , whereas no member of a different orbit may be so obtained. These orbits are the 
abstract nonisomorphic structures. 


for case II where is prescribed, 

Pm ‘ ” and r *\ 

and for case 111 where ^ - c? is prescribed. 

7 


In the limit where M 


stability criterion for all cases 


Group theoretic reasoning yields the following expression for the number of orbits or structures of a 
general m-adlc relation (the matrix becomes an m-dimensional array rather than 2-dimensional) 


» ,r * s? X. fi 7r) 


where str K is the number of structures of m-adic relations on a set of n elements. To define f (rr), look at the 
result of a given permutation TP upon ea h of the 2 n relational arrays. Some of these arrays, as pointed out 
above are invariant under TP, that is 


In the case of heat transfer in a finite slab O * X i J when temperature is prescribed on the boundary, 
the explicit formula of finite grid network is stable for r values greater than - depending on the grid size, 
but when a heat transfer f.lm coefficient is present, the formula is stable only 2 for r values less than * . 

This latter result cho« ks with Fowler [2] who derived the similar result by contour integration method. 2 


ft ft*) is defined as the number of A's for which this holds true. The summation is then taken over the n! different 
XT" s. On the surface this formula would appear to involve inspection of all the n* transforms of the 2 nm arrays, 
comparing the n m elements of the transforms with the elements of the original arrays, hardly a feasible project. 
However Davis develops a version of the formula applicable to computation: 

m 1 


Andrew T. Ling 


References: 


[ij W. E. Milne. "Numeriial Solution of Differential Equations", John Wiley. N.Y., 1953, p. 122. 

[2] C. M. Fowler. "Analysis of Numerical Solutions of Transient Heat-Flow Problems". Quart Appl 

Math. 3. 361 -376 (1946) 

CALCULATION OF NUMBERS OF STRUCTURES OF RELATIONS ON FINITE SETS 

In group theoretic studies, the question has arisen, given a set of objects how many structures of 
relationships may be defined among members of the set? A general formula which yields. In principle 
the answer to this problem is available, but even for small sets this formula involves extended calculations, 
practically demanding a htgh-apeed machine for their accomplishment. Whirlwind 1 is at present being used 
to evaluate the formula for small sets. 

To specify the problem precisely, consider a set of n objects. A dyadic relation among members of the 
set is specified completely by an nxn matrix of 1> s and 0 - s. a one in the ij place indicating that element i bears 
the relationship to element J while a aero indicates the absence of such a relationship. Counting the number of 
structures of relations amounts simply to counting the admissable arrays of 1' s and O' s in such a matrix. 

With no further restrictions, we see that the answer is 2» 2 , but in this figure of 2n 2 we have included many 
isomorph.c structures which can be permuted into one another simply by renumbering the objects of the set. 
iru» then va the problem, how many noniBomcrphic structures exist? 


, / ,r, . f 

b (ir) - n ! ( 1 f> 1 l 


or in particular 


c * least common multiple of r^. 

(h, k) * greatest common divisor of h, k 

There is still plenty of calculating involved in these formulas. Consideration of different Tt 1 s has been reduced 
from all of the n! permutations to representatives, 7^, from each of the conjugate classes of permutations 
having similar disjoint cycle schemes 


Where P r is the number of cycles of length r in the permutation. (In a cycle, the first element is permuted into 
the second, the second into the third, the last into the first). Since the total number of objects is n, the must 
fulfill * 

X> Pr 

V • I 


i 
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Mont of the p are aero. The total number of conjugate claaaea. or terma in the summation for str^. I* the 
number of partitions of n into Integral summands. For n • 10 there arc 42 terms for n • 20, 385. 

(Compare these figures with nf, the total number of permutations which is 3.6 a 10® for n • 10 and 
2.4 a 10*® for n • 20.) The redundancy, or number of member permutations, of each conjugate class is b (*r). 
f ( tt) has beer, replaced by i'-^where d represents the number of ’’degrees of freedom" corresponding 
to a particular permutation. The number oF degrees of freedom is the number of elements of the array A which 
one may choose at random and still be able to have TJJiA » A. This number is definitely limited since in 
choosing, say, ojj we must assign the same value to Affine* and in turn to r*< tr^.etc. until we go through 
an entire cycle returning to s^j. In the particular case of dyadic relations, m • 2, the computing formula for 
djCfr) includes a term for degrees of freedom in the off-diagonal elements and one for the degrees of freedom 
in the diagonal elements. Other formulas are available for special classes of dyadic relations, symmetric, 
reflexive, or antisymmetric. 

This problem is interesting from the programming standpoint in that it requires exact answers which 
are very large integers, far exceeding the 15-bit capacity of Whirlwind I. Operations must be conducted in 
integers with no round oft; in division, remainders must be kept. The program is arranged to develop each of 
the position schemes (p^, . . . p^l in logical order. From each partition follows b (7H) and d^(7P). Extensive 
use of the logical orders is necessary in the partition developments, and in efficient calculation of the sums for 

d ITT) which, as written contains about n^/2 terma for on nxn matrix, but at most about n non-zero terms, 
m 

Because of the necessity of carrying very large numbers exactly and because upper bounds for the 
quantities involved vary by as much as the capacity of several whirlwind registers even within the calculations 
for one particular n, (for example upper and lower limits on b (tt) one n! /(n-1) and 1) and are in any case 
difficult to estimate, a flexible arithmetic routine has been provided which is capable of carrying out addition, 
multiplication, division, and shifting of multi-register positive numbers of arbitrary length. Opciation with 
this routine will be very slow for two-register numbers, but unfortunately numbers as long as 8 registers may 
be encountered even for n »10. 

The present calculations are limited to the calculation of the numbers of structures at dyadic relations, 
and will be carried to sets as la»-ge as time will permit; 1 hope to n *= 15. The necessary routines are now being 
tested on Whirlwind I. 

Reference: [l]R. L. Davis. Proc. Am. Math. Soc. 4. 486(1953) 


M. Douglas Mcllroy 


THE STABILITY OF THIN, SHALLOW ELASTIC SHELLS 


The problem under consideration is that of the stability of a hyperbolic paraboloidal shell loaded by its 
own weight. Specifically we wish to find that load, p which causes the shell to buckle. The equation of the shell 
is ° 


(1) / - C * i 

a b 

and the assumption of shallowness is that ( ~ ^)*« 1 • 

The precise conditions of the problem being studied are: 

a) The shell is uniformly loaded by p . 

b) The edges of the shell at x • 0. a and y » 0, b are assumed to have moment free support. 

c) The edge stiffeners of the shell are assumed to be rigid in the direction of their axes and to 

have negligible bending resistance in planes tangent to the shell. 

Using these conditions a complete solution of the general equations of shallow shell theory is il] 

(2) w • 0 F • - 

2c ^ 

wher*. w l- deflection In the a Jiievtum mad r is the stress function. To test the stability of this solution we 
assume new solutions of the form 

12 
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(3) w • 0 ♦ w 

Substitute these into tnc diffe* 
The equations which held are 


r * ' B ic ♦ # 

•al equations and linearize the differential equations in terms of w and 4 > . 


(4) 


V l v l f 

D v'v * & - -If f *4- t- t> -4. 

where h is the uniform thickness of the shell, E is the Modulus of Elasticity and D is the modulus of Flexural 
Rigidity. The boundary conditions are 


(5) 


x * 0 J a 
H ~e>, k 


^ - V*£> 


- 


(4) and (5) then represent a characteristic value problem for p Q , the buckling load. 

To find the buckling load we assume series expansions for w and . satisfying the boundary conditions, 
as follows 


2Z Z *• '-P ** sj* 


( 6 ) 


svx 




*.•< *«• 

Substitution of (6) Into (4) leads to a doubly infinite set of linrar homogenous equations in the A mn and B . 
P is a parameter in this doubly infinite set. We desire to find the smallest value of p Q which permits solution 
of the set by non-zero A and B . To obtain approximation* to the desired p the series (6) will be terminated 
after a finite number of^erms. solution of this finite eigenvalue problem and the convergence of the solution 


to the true value of p q will be investigated. 


References: 


Anthony Ralston 


[l] E. Rcissner. "Cm Some Aspects of the Theory of Thin Elastic Shells." Journal of the Boston Society of 
Civil Engineers ( I* P*est) 

MACHINE SOLUTION OF THE DIFFUSION EQUATION 

The one dimensional diffusion (or heat) equation of the type 

“» If * • kc "' 9f "■ °- x * L 

with initial and boundary conditions 

-k(x) §£ | *&,+*■* 


* • O 

2 JL I 

d X ' M • L 


o , 


U(X,o) 3 0 , 


o < X < L 

where U is concentration (or temperature) (X has units of length and t has units of time) has a unique solution 
which can be put into the form 


(2) 
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wtar * * I ( [^ 

^’ -? fzj*> ( /^; f zj*) s. . 

9 O O 

The Xh are a complete set of orthogonal functions belonging to the associated Sturm-Llouvllle system: 

£i k( ’> 7 # k J * =* 

u r 

^ * in.* ^ * I x.£ 

This solution is found by using the separation of variables method. First, assume the form 

U* 1 * R (X) i S (t) ♦ C and satisfy the first two boundary conditions of (1). This will serve to determine 
a constants except This is the steady-state solution (that part of the solution which does not die out with 

time). Next use U » X (x) T* (t) with the boundary conditions for (3) above plus the condition on T* : 

| K *r */f) _ $ This product form yields the solutions 

t 1 * -X.t 

c k e x/*). 


These two solutions are then fitted together (using a type of Fourier expansion of R(x) in terms of Xn function) 
to satisfy the initial condition at t ■ 0. The constant c is determined from the fact that 


Im 

flt(‘, t) Jt * Qt. 


(See Quarterly Progress Report No. 14). 

Equation (2) is valuable for two reasons. First, it displays the steady state solution and second, it 
suggests when the transient has become insignificant. A simple approximate calculation showed this occurred 
when t « O (L). 

fl o 1 

A number of cases have been solved analytically in the literature L ' J , When K is constant, Xn are 
sines and cosines; when K is linear in x, Xn are given in terms of modified Bessel Functions of zero order. 
Corresponding forms of K will produce Legendre polynomials, Legarre polynomials, etc. 

An interesting case is when K(O) r 0. For instance K (x) - x. Thus, by the Laplace transform with 
respec t to t equation 1 becomes 

x ** ** V i* 

7w * t— - S V • 0 
j x jx 


— f - it 

u(*,s) ■ J e Jr 

o 

The solution to this is readily found to We 


^ ^ (2 y/sj, ) where is the modified Bessel Function of zero 


order. The inverse transform then gives 


j^c sr - - a fi(- */*) 


Two cases of K(x) were solved on the IBM Card Programmed Calculator. The first was for a constant K. 
The cpc solved the transient region in a little over two hours using the difference equation lc>: 

U (x, t + k) - e U (x + h. t) ♦ (1 - 2 f ) U (x. t) » > U <x-h,t) 

Kk 

where r * -j^r . k * t spacing, h ■ x spacing. The h and k were started very small, then, as the calculation 

proceeded, h was doubled and k increased by a factor of 4 so as to keep the same difference equation throughout. 
The h was doubled in this manner 8 times. The final solution agreed with the steady state solution to within 0.1%. 

The second case was solved using a K as shown in figure 1. 


K (O) << K (L) 


K (-b) - O 


where b s.< l 


In gene-al, the solution U has a singularity at X * - b which causes considerable difficulty in approximating the 
boundary condition at X * O, since any power series expansion on or near the boundary necessarily has a small 
radius of convergence. This seemed to suggest a variable spacing in the x direction (and consequently in the t 
direction). The difference equation corresponding to the point pattern in figure 2 is 



(4> - 4 *t r C. 

where A„ - (.K.' )/' ^ ’ L/V ^ /C,' ) . \ r L, ) 

ttnd T3 h = 1 - -!>..,)]/k H A._. 

where K„ - K(*J, X k . jt A, ; I = k'C*)\ ^ 

t 9 * * V - X* ■* w \ 

It will be noticed that An + Bn + Cu « 1 and for local stability, all coefficients should be positive, consequently are 
equal to or less than one. Thus, we have the two conditions: 

(s) 

For constant K and constant sparing these reduce to the well known stability conditions 

M. A.‘ - kk/h* < '/t , k > e>. 
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The first one In (5) places the limitation on h n while the second, for given h n and h n ( gives the restriction on k^ 
It was found that the x spacing satisfying the relation 


(6) 


f- C i *. i L 4 (*> l) 


will conform to the first condition in (5) very well. Certainly, such a pattern is not unique and would depend on 
problem as well as the inclinations of the individual. This particular one worked very well for the KQ that was 
used. For the rest at the interval (1/2)L 4 X n 4 L constant spacing was used. Thus, for the first half 
interval 

h • 2 h . 
n n - 1 

The kn spacing was chosen such that k -2k . for the first half of the x -interval and constant k was used 

^ n n -1 

thereafter. i.e., for *n > ( ^2) L. The following figure 3 gives an abbreviated version of the overall point 
pattern used. 


-p I 


: / 
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Even with this short-cut of variable spacing the calculation time is very large due to the relatively dense 
columns of points for n • 0, 1, 2. This is especially true when relation (6) holds for a dozen or so points before 
. onatant spacing is applied as was the case actually used and described by figure 1. A further reduction in time 
may be had by combining the difference equations <7> in an appropriate manner. Thus the first equation with the 
above boundary condition at X • O is 

* ^ a A(.) * (*.+TS>) - C K ~ 

and U, t « \a/\ct») * (e X) \ t Q *..« 

= k. a/kM + /w«i . + - c k,J] 

♦ c. 1-4, V - C.u.'_l t 

or W. e a, * (. u ♦ J u, «. e if. _ 

“ ^ * 9 0 * h • / • j 

* i . * k - r 


likewise U, 


,) _&a_ 

: 7 * ’•* • s 


A •" 

/ 




where 

a 

o 

■ Tw> <l * A o* B o’ 

•i 

. a 

I K(0» 



>4, 

r ! 




b 

o 

■ <A 0 + B / ♦ A, C o 

b i 

* A ,<VV + B i 

A i 


J h* --J 

_ hilAx _ 

! _ 


d 

o 

■ (A ■* B ) C ♦ C B, 
o o o o 1 

d . 

* A . C o * B l' C l 




The dotted line represents the maximum extern one may predict into the region from the points on the bottom 
row (solid dots). The points calculated from the row are shown in open dots. The rest of the points are calcu¬ 
lated using more immediate points (the open dots and previous crosses) and a completed line of points is eventually 
obtained across the entire interval (0,L). 

The difference equations would be as follows 


C. C 
1 o 


One may proceed in this manner for U 


B . c , 


C 1 B 2 


C C 
1 2 


U 


u. 


u. 


u 


o, m + 4, o, m ♦8, w l, M t- 4, “1, m ♦ 8. w 2, m ♦ 8, , , . , as long 
as one likes. The coefficients will always be in terms of sums of products of preceding coefficients and the 
A. B, C. The resulting point patter may look like that In figure 4. 
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where , . 

A. *./ t C. * C.U T5.' , 

The boundary condition at X • O may be approximated by 

<« "w- « - *.«/>Cfej. 

The condition at X = L may be approximated by 


U. 


N-l.m 


N+l, m. 


where X 


N - 


ht t r —* 

A good checking device may be had when calculating these coefficients if it will be noticed that all 
coefficients for each U must sum to unity (with the exception of the constant term which may be present). 

For a dozen or so sets 6fsubStitutions yielding + 40 98 U + <096 etc., one may resort to a high 

speed computer. It was found in this problem (the one described by figure f) that combining the difference 
equations in this way and using the CPC to compute the coefficients, the total overall computation time was 
reduced by 20 fold. 

A word or two is in order concerning the boundary condition near the singularity. The solution U changes 
very rapidly in this neighborhood especially for the transient solution and the boundary condition (8) is not 
completely adequate. It was found in an early run on the CPC that the approximation (8) introduces a very nearly 
constant negative error into the interior of the X f t - plane. This had the effect of using a different reduced Q 
so that the condition: 
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was not fulfilled. The ratio was conaiatantly leaa than unity but essentially constant. The boundary condition (8) 
was used in calculating the coefficients, then when the problem itself was ready to run the boundary value ^ 
was calculated from the condition 

L 

- a* • • *!**•,*» ♦•••*• M 

b 

i.*.. the relation was solved for U . The M ' s »re the weighting factors for the numerical Integration. This 
la essentially equivalent to the differential form of the boundury conditions so It la a valid approximation. The 
error Introduced by the numerical integration was smaller than that in the difference form of the boundary 
conditions (8) at X • O. Also, in this integral form the errors tended to be self correcting as the calculation 
proceeded. The U ,U^ ^ values obtained in later runs were much higher than one would expect, but all of 
the real of the poinVtA’ • quickly aettlcd down to within 1% or leBa of their true values as was shown when the 
steady state solution wife reached. Even this error could be reduced by taking the mesh in the x direction finer. 

This method of combining the difference equations is subject to considerable generality and appUcability 
to other problems. It combines the higher accuracy of very fine mesh spacings (provided finer meshes provide 
higher accuracy to begin with) with shorter solution time found in more coarse meshes. The number of numeri¬ 
cal operations is reduced by as much as three or four orders of magnitude for very large problems i.e.. ones 
with difficult K(X)'s to handle, thus reducing the ever present danger of accumulated round-off error. The major 
limitation to this procedure is that of calculating the coefficients. 

When the solution time of the given problem has been reduced to the point where it is roughly equal to the 
time required for the computer itself to generate the coefficients, It U felt that this is probably the point of 
maximum efficiency for the method. 


Philip L. Phipps 


References: 

[ l] Carslaw. H. S., and J. C. Jaeger: Conduction of Heat in Solids, Oxford, The Clarendon Press. 1987. 
[2] Churchill, R. V., Modern Operational Mathematics in Engineering. McGraw-Hill Book Co.. 1948. 


MULTI-CENTER INTEGRALS 

The problem of determining multicenter integrals of atomic orbitals Is still fairly difficult. The following 
report will review briefly the status of numerical procedures with respect to the Whirlwind computer here at 
M.l.T. For general references the reader should consult a recent review article by Dalganci (1) . 

One method of evaluating molecular integrals which is specific to two-center integrals is to limit the 
atomic orbitals to Slater AO's (Le. analytic orbitals of the form |i) v" 'e'**’ 5 (6 <f>) 

and then by use of prolate spheroidal coordinates express the integrations in terms of analytic functions. When 
this is done it is found that all one electron integrals can be simply found in terms of the two functions A (X) and 
B n (X). which are only sparsely tabulated. The latter functions, however, can be generated for arbitrary arguments 
by means of a high-speed computer and such a generation subroutine has been written and tested[2L (In addition 
tor convenience, a production-form program has been prepared for all onr-clcctron Integrals of orbitals up to 2p; 
computation time about 1.75 seconds per Integral.) The remaining two elcctron integrals fall into two groups, 
those with the two charge distributions each about one of the centers, (i.e. coulomb-like), and those with split 
. barge distributions, (I.e.. exchange like). The coulomb-like integrals are also expreasable in terms of the A 
and B n function* and thus offer no great difficulty. The exchange-Ukc Integrals, however, must be evaluated n 
by different techniques, which have been deacrlhcd in detail by Ruedenbergpjof Mulllken' a Group In Chicago. A 

Gf*nunhand )° e ” Iu ‘ ,le th<, “ ••change like integrals has been written by P. Merryman. also of Mulllken' a 

G. oup and 1. currently being tested. Thus all the two-center integral, will soon be fairly easy to compute. 
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There are two major limitation*, though, to the method juat described. The first i* the requirement that 
the orb itala be Slater AO' a. For Hartree - Foch orbital*, thia problem can be circumvented by fitting the function 
by mean* of a few Slater AO' a [4.2]. Beside* the difficulty of fitting, the only other disadvantage of this decompo¬ 
sition la the greatly increaaed number of Integrals required and the consequent many numerical recombinations 
necessary to obtain the basic integral values. The second limitation of the above approach la the obvious inability 
to handle three-and four-center Integral*. 

A aecond method of calculating multicenter integral* exist* which does not have these two drawbacks. 

This la the method of expanding all the orbitals about one of the center*. Thus If the orbital* being used are 
expressed in the form, f (r^) S/* ( <?*) about center a. they can always be expanded about another center b, 

by virtue of the completeness of the spherical harmonics. 

f(ftV 

(In addition. - and £ can be expanded In the usual series of spherical harmonics and powers of r.) Now 
when f (r B ) i{ & numerically-given function, as Lowdln has shown [5j, one can use straight-forward numerical 
procedures to find the expansion functions, g (r a p, r*,). When the basic orbitals are expressed in terms of 
Slater AO's. though, the function g la even simpler to obtain, for In this case, the g's can be written In terms 
of the spherical Bessel functions of Imaginary argument i (x) and k (x). New numerical techniques-have been 
found for easily generating these functions on a high-speed compute 1 }*, and such a generation subroutine has been 
prepared for Whirlwind. In any case, the expansion functions, g. may be systematically prepared. Furthermore, 
since the multi-center Integral has now been reduced to a series of one-center integrals, the integration over the 
angular coordinates can all be performed, a process which greatly reduces the various summations introduced 
by the expansions. Thus the basic integral is expressed in a series of terms each containing a radial Integral 
which will be either of the one or two electron form: 

00 

o 

<*■ [yfen iio^vA 

where a, b, c, and d are, in general, numerically given functions about one of the centers. Because, in a complete 
calculation one may require several hundred of these integrals, a Whirlwind program for doing these numerical 
quadratures has been written and tested. In the interest of program Bpecd, Simpson's rule is used for the 
integrations, but the mesh size can be essentially arbitrarily small and with as many scale doublings as desired. 
To give a very rough notion of the computation time required, a basic 2-outer exchange integral might involve a 
series of about six significant terms, each of which requires about five seconds to compute, giving a total time 
which compares well with even the explicit two-center method. 

Thus the method of expanding all functions about one-center appears very attractive for the three and 
four-center integrals. The writer plans to use this method to evaluate the more important three-center integrals 
occurring in a tight-binding calculation of Graphite. 


F. J. Corbattf 
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EIGEN VALUES FOR A SPHEROIDAL SQUARE WELL 

L For conciseness the description of the problem will be restricted to the case of the oblate spheroid [i]. 

We define the co-ordinates //, which are defined In the regions / * S' * I The problem is to 

[ o* * /rr J 
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, A >* ■ 


*» <+ 

The ^ and are conatanla (to be determined) and the 5 h C mi functions are as defined in Morse 

and Peshbach. m * *9 j 

In order that the V which Is of the form (2a) for u< <L and of form (2b) for s^bea solution of (1) at 
* 4 . It Is necessary that 

, 3) fV. 

K^.4. ‘ ° 

Defining the function of g and h 


wtar /*' " rA w f t 

I S *i ( p' s *.„ < /, » J' Ll/ , 

the equations (3) lead to the matrix equation 


(O A • a • O 

where ill! vector with component* «, and A la an Infinite square matrix with element* 

In the usual way it U seen that a non-trivial .olutton of equation (4) may be found by requiring that 
det A • O 

and the values of k* . h / ^ for which thi. determlnental equation la true are the eigen value* of equation (1). 
If one choose, to think of the C,, In term, of their role a. expansion coefficient., l.e. 

■ 5 -' Kh,>> * X Xl C L* S K X ( Jj>) 

it I. immediately obvious that in general the otpan.lon on the right cnnot be expected to terminate For if it 
relation.nip between ^ ( kjJ and the p y . But such a relation.hlp doe. not exirtjl] 

On the other hand the differential equation for \ t which 1. 


I X ~ - Ay/ 5 *, -■ 


XXIXXX >M ttat for A.,» the differentia. 

F ' or <orvlct.v.r..),heX ? ^rX.dr.^Xr,h.f ^ 9*-° 


11. It seems likely that it will be necessary to compute a large number of spheroidal functions in order to 
solve the problem outlined in the previous section. In order to do this 1 plan to use the Whi-lwind program pre¬ 
viously developed and reported upon by F. Corbato and J. D. C. Little.[3] However, Mr. Corbato has found 
that his program in its original form could not be used for values of h much in excess of h>8. For values of 
h much in excess of this value the iterative process used in the original program becomes extremely sensi¬ 
tive to the initial guess (for the eigen value), and if this guess is not sufficiently accurate the process will 
converge to the wrong eigen value. In order to improve the initial guess we (Mr. Corbato and myself) are 
computing the eigen value directly by means of an asymptotic formula.f4] 

Preliminary runs seem to indicate that the asymptotic formula gives a surprisingly good estimate of the 
eigen value for values of h as low as 5. If this behavior is consistent it may reduce the average time needed 
to compute spheroidal wave functions by as much as 30%. 

Jack L. Uretsky 
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VARIATIONAL DETERMINATION OF ATOMIC WAVE FUNCTIONS 

During the past quarter, most of the time has been spent programming (or Whirlwind. The variational 
calculations for the three parameter wave functions of Morse, Young, and Haurwitz. [1] 

The actual evaluation of the energy for any particular set of parameters is very easy to program. There 
arc two basic integrals, the first being the rather trivial function 

<0 A h (*) * f - i'.'/*'"*' 

O 

involved In the matrix elements for the kinetic and nuclear potential energies. The second one is involved In the 
electronic repulsion and exchange terms. 
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( 2 ) 


f>i r C~P r ' r ' (c*. <?,) 
It”* f (hi-1*1)1 ] 

U» * if l (*- 1 * 1 )! J 


A yjc f ' Ic.s <fj) 

sjfc f 


I 3 

The function S ^ is in Quarterly Progress Report No. 9. September 15. 1953. In general, 

it can be written uh 


Sl(f,r.*) * (t *+**'*'■) /**(**$)*' , * - ■ 

K* 4 ? 

and are integers and along with the » r are functions of p. q. and )i . 

The central problem has been to obtain, using the shortest possible amount of machine time, the energy 
minimum correct to three decimal places. 


The method has been tried for the ground Htatc (/$* <IS * P) * V of Boron, and Ihe results appear to 
be correct to one unit in the third decimal place. Several further checks are being made by starting at different 
points and using different values for A^ , If the results are satisfactory, a systematic program for calcula 
ting the energies of the states not yet computed in previous hand and IBM calculations will be set up. 

The possibility of adding a 3s function to the present scheme in mow being investigated by two seniors 
Robert Papa and Alfred Finn who are doing a joint B.S. thesis in the Physics Department on the calculation 
of the ground state of sodium and several other states involving a 35 orbital. They are working with a three 
parameter form for the 3 S function with the hope that they might deduce .i rule fur fixing one or two of the 
parameters. 

Finally, programs for calculating the configuration interaction integrals Involved in the pertubation 
corret tion terms‘^Iare being written for substitutions involving JTS, ip , and 14 orbitals. 


Arnold Tubis 
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The first minimization scheme tried was a second order Lagrangian interpolation. It proved to be highly 
inefficient, a large number of points over a small region being required for a reasonably accurate extremum. 


NEUTRON-DEUTERON SCATTERING 


During the past three weeks, a new approach has been tried. The technique, commonly referred to as the 
method of steepest descent, may be summarized as follows. An initial guess of the parameters (a . b , c ) is 
made and the valve of the energy and the partial derivatives & W are determined °From°the 

d A ' ? L J jo. 

latter we obtain the gradient vector whose negative direction is along the path for which the energy decreases 
most rapidly. Since we are looking for a minimum, we take the next point along a line parallel to the gradient 
vector. By properly sraling the distances between successive points, we should be able to converge to the 
minimum. 


For the actual machine computation, we first estimate Hie maximum distance of our initial point from 
the minimum. Let this be called A^. The next point is then chosen to be 


The variational calculation of S-state phase shifts in neutron-deuteron st uttering is nearing completion. 
We indicate here an outline of the theory. 

Assuming two-body central forces of arbitrary exchange Ivaracter. the Sohrodinger equation for a 
system of three particles of equal mass takes the form 

where 


Co: cPj 


+ V 
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h ! * b Q - C«5 <y b 
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U i j <* + n> 




♦ b B, 


I. M 


V 


e l • c o ■ 4^ ^ «fc 

where Co, .etc. are the direction cosines of -he gradient. After that four more points are calculated 
the distance between successive pomts being acaled by a factor of two each time. Such a procedure should 
give an extreme point with the maximum uncertainties in the parameters given by 



(C.s /it 


A, (fescPbL.» / It 


A 0 U-i <?c / u 


where the direction cosines are evaluated at the point of minimum calculated value. 

The partial derivative, are calculated by taking first differences of W for parameter Increment, of o.ol. 

u*log the generalized ronr.pt- of distance, gradient, etc. In N dimensional cartesian space we may 
* r ,d the preceding argument to the caac of functions of any number of parameters. 


« and Btj being the operators which exchange the space and spin coordinates, respectively, of particles i 
and j. Ui} is a scalar function of the distance y and w, m, b. h are numerical coefficients whose sum is unity. 

After the motion of the center of mass is separated off, there remain six coordinates tojlcscrlbe the 
system. For the scattering problem in question the appropriate choice is that of the vectors >* and y , 
where P Is proportional to the deuteron coordinate and f ia the vector from the Incident neutron to the center of 
mass of the deuteron. (Fig. 1; it is assumed that particles 2 and 3 form the original deuteron pair), it can be 
seen that the magnitudes r and . together with at, the angle between the two vectors, are sufficient to 
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determine the relative position of the three particles In their plane. The remaining three coordinates then 
describe the orientation of thia plane. For low energy (S-state) scattering, the wave function will be independent 
of thia orientation and will thus be a function only of h . / and a . These coordinates are of course expressible 
in terms of the uiterparticle distances, which furnish an analogous specification of the relative location of three 
particles in their plane. The relations are: 


if y 
l r *> 


13a) V- 

( 3 b) = Vi * ft 1 ) ~ 1 j// 


( 3 c) 


f (ft 1 * ft 1 ) - t/} . 


Instead of ol . we shall eventually use as a third coordinate the distance S, the third side of the triangle formed by 
and . This is given by 

(3d» s - + 

We also find it necessary to define the coordinates y , y the result of /'•^j (neutron exchange) operating 
on the original and . (Fig. 2) \r* and s' are expressible as quadratic forms in * * J and S. 

,4 " > K ' m })T(f*(3-fj) r r^t-zr) t Vss*) 

Mb) -f' -- j **(1*ST) -as*) 

Me) 5' * $S‘) 


Writing equation (1) in the V*, ^ system and expanding in partial wbvcb, we have for the S-wavc part of 
the resulting equation: 

where E K la the kinetic energy of the incident neutron In the center of mats system ( • | the energy in the lab 
system) and E^ Is the binding energy of the deuteron. ^ 

The wave fun. lion if in (5) is a combination of spatial and spin functions, and must be antisymmetric in 
the two neutrons to satisfy the Pauli principle. 1 We therefore write explicitly, for the quartet and doublet 
spin states (S « 3/2. 1/2 respectively!: 

<6.1 y a - Xa (id - ) 


We do not utilize the isotopic spin formalism, ss we find II more convenient to label each particle expllclty i 
neutron or proton. * v 
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Here \fu is the (totally symmetric) quartet spin function; and is therefore multiplied by s spatial fun. lion anti 
symmetric in the neutrons ( ^ In the doublet state there are two independent spin function* (for each 

value of >. We choose J(p j to be the one symmetric in the neutrons and fti the one antisymmetric, and 
multiply each by the appropriate space function so ss to sntisymmetrixe the total wave function. We now put (6a) 
and (6b) into (&), introduce the notation y ■ jfp end sum over spin coordinates. The resulting equation for the 


quartet scattering is: 


[ T f j|_ ~ kj ](</-*) + rq [u -u )/y, 
y i J 


( 7 ) L 1 - i («-*} + >-?-U lit-*'Sty *£> 

where we have introduced the following additional notation: 


( 8 ) 

(8) U 


T * 'rj t rP * ?f' * ( * C0t * 


*T " 


,10.1 f. tr lfl 

uob, It/ = fre, 

doc) Cl - - V 

7** 

For the doublet scattering we obtain instead of (7) the two equations 

( 1 U) [T - k ~ jt ] ( « - S) + Y f[ i- j * 

u yp j > v<? ** yj 

[ T- h!^‘] (u + s) w n* - 


(lib) 

where 


(12a) Ojj c U u [ u b * (t*+k)to, t ] + i'Jj [ (<0 - £■) + (*• ~£r) to,i ] 

* Vu * <*•-£) ■ 

<I2W Jia - ^W-1/ If £ 

“ic* ^kb - V u [uj- b * (» -11/4,1 J + *£)*("*£) to,s 

We shall proceed with the derivation of the quartet scattering. The doublet gives longer equations, but proceeds 
analogously. 

The boundary conditions on Ip and ^ are the following: 

a) J must be everywhere finite; hence u • & when either Y ■ O or ^ ■ & . 

(13) b) when f+*o, V -+ W (r) Sf*( let+S)/f ^ 
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when F (►) la ths deuteron (round state wive function, and 5 in the phaie ihlft which deicrlbea the scattering. 
Hence we have, for J -+ oo y Y < oo t 

(14) U >Fl * ) St-O-f+f )/»‘- J> SW- W, -r t^cSJ 1 

(15a) where Wj - >-FYV) &r(kj) 

05b) Wg • ► 

W ( and Wj. and hence W. the aayntptotlc form of « . satisfy the equation 


r r ' k Tr']l*'l + 

W{l l ^I 


± -0 


Together with (island (14). we have the boundary condition for exchange scattering, y ’ 00 J I" ^ 

(17a) y —F U* (kj'+S)/j'v* J 

- , - _ . .. n 1 / Z 1 — .-7 ■ w'. Wi faf 


(17b) « 


>-TrCt-9 - w - w, - 

r 


where now F( > ) la the deuteron function In the exchange coordinates, and W [ and W 2 are given by 
(18a) W, Lt F(y') 

dab) w t * F(r’) S>“ (kj ) 

W. W t and are solutions of 

(19a) f W ) f '*' ) 

09b) [ r -• u -j-L jJ ]{ V?, + >?■&„ 1 * t \-o 

(19c) *T [ Wj J ( #2 J 

The actual boundary conditions are of course on the total wave function of the system, which Is V 1 - V*'. 
However It la impossible for both if and to contribute at Infinity, since they go to deuteron functions in 
different coordinates ( V and >•'). That la. If one neutron la at Infinity then either V or V- must be Infinite and 
the deuteron function, which drops off exponentially will lull the contribution of that part. This means that (P 
Is used to satisfy the boundary condition for ordinary scattering and Y ’(*■ a)tor exchange scattering. 

At this point we Introduce the inside wave functions ^ • f • defined by 
(20s) tj = tO - u 

(20b) us. Q - ei 

From (7), (16) .nd (19) we obtain the differential equation for tf- f / 

m> [T. tlS/jfy-f)* + 

* > 


Tt»e boundary conditions on 

( 22 a) 

If ( f * °) - 

yF (bj 

(22b) 

<J ( V = 

o 

(22c 1 

^ ( J * <*’) * 

o 

(22d) 

^ ( v- =.*>) - 

o 


Starting from equation (21) many variational principles for the phase shift £ can be constructed, as 
pointed out by Rubinowfl]for the two-body problem. 

The general procedure is as follows: from (21) (16a) and (19a) we write the differential equation for 

(23) fTa - UJ t - 

1 ^ />•] ■* Flu u;, y>-f ]} 

Multiplying this equailon by ^ ' . (16a) -(19a) by Off- M- c w J , adding and integrating over all 

apace (/U •= Y l j S -ciHl we obtain, after a partial integration of the second derivatives, 

<«> -gf fcc.f'eT) = C 4 L t 

.0 «, 7? 

o o * o ' * 

+ (k‘ - \.j l )(r f >*♦ f-f <(y--f * o,)n(p. u ^.i*3 )+ (if- S- ■?y • ^ - V. ^>} 

(25c) 4 similar derivatives in v* and ctf 

Again, on multiplying (23) by (W - W ), (16b) -< 19b) by ( <f *- ^ ), subtracting and integrating 

we obtain 

(28) B -+-L, a A k c.i S' 

where 

(27a) A- j~ h l J 7 r 1 [^ 4 -tP 4 )L (n-n-jHt/bf ~ 

(27bi Lj 1 / J* f I ^ r' ^7 ( u * ) If" 

J 6 4> > 

4 similar derivatives in V* and <* . 

Lj and are determined from the boundary conditions, and do not depend on the nature of ^ . Neither equa¬ 
tion (24) nor (28) is stationary, but suitable combinations may be taken which will have the desired property of 
stationarity. The simplest of these combinations is obtained by rewriting (24) as 

(SB) (13 + L t ) l< cSS - C + L l r L, ( t eSS) 





















APPROVED FOR PUBLIC RELEASE. CASE 06-1104. 




GRADUATE SCHOOL RESEARCH 

, ;•> 

Substituting (or k'.uTV fruni (26) we obtain 

(291 h , J J - ! [ ( B r A,) A - ( C* i-,) ] 

L l 

It may be verified that (29) is stationary with respect to variations in «• and which satisfy the boundary 
onditions (22). The problem now consists of choosing a suitable trial function J . dependent on a number of 
variational parameters, and "minimizing" with respei t to the parameters. Since we are dealing with phases, 
which are not positive definite, it Is not known a priori whether the stationary point will be a maximum, minimum, 
or a point of inflection. 

The integrals which appear Ui expressions <25) and (27) are in general extremely complicated, mainly as a 
re mit of the exchange operators and the mixture of primed and unprimed coordinates. Any function linear in one 
net of coordinates is a sum of radicals In the other set. and the exchange operators introduce various sets of such 
radit als. The best chance of obtaining a reasonable solution lies in the use of Gaussian forms for both the 
potentials ami the trial functions, since linear forms in the squares of one set of coordinates will again be linear 
in the other set. Accordingly we > house for the spatial dependence of the two-body potentials. 


where we have taken the range of the force as our unit of length. (All the calculations are performed in dimen 
sionless form). The depth V Q is a compromise between the values required to find the deuteron and the triton 
■ orrectly, in the absence of tensor forces especially, these two data are at variance: a well deep enough to give 
the ■ err* * t binding energy of the deuteron gives too much binding for the triton. The values of the constants 
w. m, b, h in (2) are left arbitrary, except that the singlet-triplet ratio in the 2 body S state gives a constraining 
relation 

(31) U/ - m - j, - h J-aiJ 3S 
For our trial function we take 

(32) ^ S '1 

which satisfies the boundary < onditions provided X + - is a constant determined by the deuteron wave function. 

We therefore have essentially a two-parameter trial function . Putting/** o corresponds to the no polarization 
assumption, t.c. assuming that the wave function is everywhere a deuteron function in one pair of coordinates 
multiplied only by a fun. lion of g. This assumption has been made in all previous calculations of the problem./2.3,4] 
The size of the "best" value of U . as well as the effect on the actual phase shifts, will be a measure of the 
importance of polaruation in the low-energy scattering. 

2 

If we expand equation (29) In powers of k . then the zero order term will be the negative reciprocal 
a. uttering length, as usually defined, and the coefficient or k2 will be half the effective range. These are quanti¬ 
ties which can be compared with experiment. The expansion is effected by performing an analogous expansion 
for each of the quantities A, B, C. L,, L,, appearing in (24) and (26). 

With the trial (unction (32). all the integrals are tractable, using methods which will be discussed in a 
future report. The integrations are first transformed to (r, q. s) space In which the volume element ts 
rqsdrdqd* and in which the limits on the s integration go from [ rqj to r -q. The integrals which give the 
most dlfflt ulty arc these arising from the "crossed" derivatives in (25b), e.g. iS. .A small number of 

dj 2 J 

these ran only be reduced to one dimensional integrals which must then he performed numerically. The remainder 
ore all integrated exactly. 

!h. ttlculation of the quartet scattering length has been programmed for the Whirlwind Computer, using 
CS 11. and trouble shooting has been successfully completed. (Results for a Inal point have been checked by a hand 
computation). The . aloulatton of the variational expression (29) for one value of the variational parameters.takes 
43 seconds of machine lime, a , onslderable amount. We must therefore be comewhat frugal in our search for 
the stationary point, and possible methods of systematic search a.e be mg considered The search is complicated 
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u, the fad, previously mentioned, that we do not know whether we are looking for a true minimum or a saddle 
("int. The plot of k versus h for ** * o (no polarization) shows a fairly sharp minimum, corresponding 

to a positive scattering length of reasonable magnitude. However, preliminary calculations show that Dm? integrals 
are extremely sensitive to the value of . and lead one to suspect that the desired stationurity ts quite likely 
to lie a saddle point. 

Leo Sartori 

.1 S. 1. Kubinow. Phys. Rev. (in Press) 

R. Buckingham and H. W. S. Massey. Pro. Roy Soc A 179, 123 (1941) 

.3 I M. Verde. Reiv Phya Acta XXII. 339 (1949) 

l 4 j R. Christian and J. S. Gammell. Phys Rev 93. 100 (1953) 

NUCLEAR CONSTITUTION 

Some headway has been made in determining the importance of the exchange potential which arises in a 
Hartrce Foch calculation for a nucleus. 

The Hartrce-Foch equations are 

(h c ' f>; + (£ / p-; v(y-y') r, (> ')->*■') Y;(o 

- (21 / vf(> '> * Cr-r') Jy ') y, <-> - A v. 

J ‘ 

V(Y y'j the internuclear potential. 

It is the third term in this equation which is awkward because it is not in operator form and it changes for 
each particle. The nucleus we are working with has 184 particles, so the solution of the HF equations would 
present quite a formidable task. 

Fortunately, the third term can be put into operator form and averaged over all particles, so each particle 
seems to be traveling in the same effective potential, or 1 pot*. The equations then take the form 

[ + L fit'*ar-)V(r-f)Vi(r')Sr' 

l 1 1 , , 

- £ ( Vlr-*') V. CVJ Vj (t •■) <jy'/?. I ifc-O ■)! I = t c O') 

4 } l 1 / 

It is obvious now that the pot is formed by the last two terms on the left, the classical the exchange 

potential. 

- » r * 

The integrations have been carried out for V ( >* - v ) * V.<? and spherical square-well eigenfunctions 
for the V ;. It has been found that the exchange potential is about 15% of the depth of the classical potential. 

The shape of the two is similar, the tail of the total pot being about the same as the tail of the internuclear 
potential. This last remark was verified by carrying out a calculation using a Yukawa interaction which has a 
long tail compared to the Gaussian interaction. 

It Is to be noted that there is nothing self-consistent about the calculations so far. The next step Is to 
find the best harmonic-oscillator wave functions by using the variational principle. The nucleus will, of course, 
collapse. To prevent collapse. Major an a forces will be added to the internucleon potential; V * V (1 ■* £ PC YY')). 
P < Yr) exchanges Y and ►*'i.e., attracts for symmetric and repels for antisymmetric states, and f is the 
strength of the exchange force, f will be varied until the binding energy of the last particle in the 1 pot' is the 
same as it was in the square well. 

M. Rotrnberg 
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COULOMB WAVE FUNCTIONS 

The irregular solution together with the derivative* of both and regular and irregular solutions has been 
successfully programmed, via the integral representation. It will be recalled that the Coulomb Wave Functions 
are given by: 

f l - 

°L r C U tyf 1 " tutor* 

V » |5 * c L <.i>{u.vf L 4 L * f Lr I 

°l’ * c L (»{ u*i) r L j L * r Lr '^'} 

the integrals being 

«o 

o 


The only essentially new integral is the one 


fjf c-l*) L c ri 


* ip f 


Everything else being done by a simple modification of the program for <$ L (y / 

The above integral 18 evaluated by a seven point Interpolation scheme, with a mesh size of 0.1. When the 

integrand becomes smaller than (10) ® the integration is cut off. In evaluating the intergrand, we "first calculate 
the quantity 

and the raiBe it to the Llh power. This procedure eliminates the possibility that the factor ( 1 •» ) 2 if 

evaluated by itself might become too large for the machine to handle. 

Some results are as follows. From our program, for the set of values (..< ~ 

<10f 2 

(10)'* 

( 10 )'* 

U0)' 1 


we get: 




F L a L-°L F L 


♦ 2.889590 
£.130964 

♦ 9.800340 
> 1.381263 
+ 9.9C9840 
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From the National Bureau of Standards Tables 

F, • ♦ 2.809859 HO)' 2 

“ -2 

F * ■ ♦ 8.131009 UO) 

G. • ♦ 9.800520 

• L .1 

O* ■ - 1.381230 (10) 

F l' °L * °l' F L ■ ♦ 1.000027 

At present we arc trying to improve the accuracy by making more use of buffers and adjusting the mesh size. 


Aaron Temkin 
Arnold Tubie 

COMPUTATION OF CHEMICAL ENGINEERING PROBLEMS OF MULTISTAGE DISTILLATION. • 
ABSORPTION. AND EXTRACTION 

The last quarterly period has been devoted primarily to interpretation and presentation of results obtained 
earlier by machine computation on Whirlwind I. No new machine programs have been developed. A few additional 
results have been obtained, largely to refine some details. 

It is concluded that this approach is very worthwhile for study of problems in this area, particularly those 
dealing with unsteady-state operations. A substantial number of qualitative and quantitative conclusions have 
been obtained. Some of these have been discussed in preceding reports. It is expected that there will be published 
in the near future. Consequently they will not be discussed further at this time. There are still a substantial 
number of unsolved problems in this field however. It is recommended that further work be done on them using 
this method of approach. 


John F. O' Donnell 

RESPONSE OF A FIVE-STORY FRAME BUILDING TO DYNAMIC LOADING 

Work has been completed on a new program which should more nearly represent the behavior of this 
building. Oir original program assumed infinitely stiff girders and thus the formation of all plastic hinges was 
limited to the columns. A hand analysis which is now available indicates that the girders in this particular type 
of structure arc relatively weak and thus permit the formation of plastic hinges. Oir new procedure allows for 
some plastic action in the girders and also computes the column resistances in each story as a function of all 
five floor deflections instead of limiting the influence to the two adjacent floor deflections. 

Our original program has been used for four loading conditions which will be checked with the modified 
program. The two programs will then be used to compute the response of the structure for about 25 more 
loadings which are now being calculated. By using both programs we hope to be able to investigate the effect 
on the structure of plastic hinge formation in the girders. A number of hand solutions will be carried out to 
check our results. 


Charles W. Johnson 

ANALYSIS OF A TWO STORY STEEL FRAME BUILDING UNDER DYNAMIC LOADING 


The upper story of the two story steel frame building is connected to the lower story by both columns 
and diagonals while the lower story is supported by columns only. We are interested in determining the 
uynamic Uwd which will cause a specified deflection ana thus result In the collapse of the structure. We alao 
want to find the effect on the values of the critical load of varying the resistance characteristics. 
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The building La approximated by a two mass eyatem aa ehown below 



FIG. I 


The upper and lower musses sre connected by a plastic resistance R I2 which prevents relative movements between 

the two masses until (R,_) is reached (see below). 
i< max. 



Thereafter « u J mM is considered to be the resistance to relative motion. The lower mass is supported by 

a spring which supplies a resistance proportional to the deflection of the lower mass but the resistance must 
t>e less than (R,,.) r 

01 max. L R o. x _ „ 4 (R ) ] 

01 1 01 01 max. 



11 


GRADUATE SCHOOL RESEARCH 


The program haa been run for six different resistanc e values and haa been checked by hand calcula- 
tlona. The data for the main part of the program la now being compiled and the problem ahould ba completed 
by the end of March. 


Charles W. Johnson 


REACTOR RUNAWAY PREVENTION 

In the past three months considerable progress has been made with calculations pertinent to the safety 
of the proposed M.l.T. Nuclear Reactor. 

To meet a deadline some undesirable approximations had to be made, but since the calculations have 
been reorganized. It is expected that the more precise answers will be obtained soon. 

During the next three months a new method will be developed which will make the calculations more 
useful and give more accurate answers. 

Little progress was made on the relaxation calculation, mentioned In the last report, due to the urgency 
of the calculations mentioned above. 


M. Truest 


2.3 Final Reports 

NEWTON ANALYSIS OF EARTH RESISTIVITY MEASUREMENTS 

Very little has been published on direct interpretation of earth electrical measurements. This has prob¬ 
ably been one of the major reasons for the periodic fluctuations of interest and activity which electrical work 
has undergone. In general, the problem of analysis of apparent resistivities to give actual resistivity distribu 
tions is a difficult one because the relations between the parameters of the distribution and the apparent 
resistivities arc not only non-algebraic, they are non-linear as well. All existing interpretation methods, to 
to the extent of this writer's knowledge contain assumptions which simplify the physical situation and, therefore, 
the algebra. The analysis presented here is no exception. Some methods also contain assumptions of linearity. 
These latter assumptions can be, at best, very poor approximations in very exceptional cases. 

The only published paper relating to direct, practical interpretation is that of Pekeris (1940), in which a 
graphical technique is developed for analysis of the Slichter kernel (Stefaneaco, 1930) under the assumptions of 
any finite number of horizontal layers, each homogeneous, isotropic, and thicker than the one above it (Rogers, 
1952). This method appears to work quite well in practical cases in which the assumptions, especially the last, 
are satisfied. It should be noted that the number of layers is yielded by this analysis. 

The method presented here assumes a given (finite) number of horizontal layers, each isotropic and 
homogeneous, with no restrictions (other than those imposed by the physics) on the values of thirknennos and 
resistivities. Mathematically, it is an iterative solution to nets of non-linear algebraic equations in several 
variables. The actual amount of arithmetic involved in any one iteration is very large, and, as many as fifty 
iterations may be required for a solution. This method would have been impracticable before the advent of the 
high-speed digital computer. Even now some would think of this no more than an exercise in algebra. However, 
only by utilizing to the fullest these rapidly increasing computing facilities will progress in geophysical analysis 
be made. 

The advantage of direct Interpretation, as opposed to indirect approaches such as visual curve matching 
or model study, is seldom fully appreciated. It stems from the excessive sensitivity of the solutions of the kind of 
boundary value problem considered here to small changes in the boundary conditions. Stated inversely, this is 
a very poor sensitivity of the boundary values (potentials, apparent resistivities) to changes wlthtn 
This can be considered as low resolving power and »• brought out from one aspect ui a discussion by 
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Morae <1993. p. 889). Thua. a very .mall difference between curve a of. aay, apparent realatlvlty can mean veiy 
great differencea In actual conductivity dlatrtbutlona. Hence. If a quantitative interpretation with any degree of 
accuracy la deaired. aa in the caae of detailed aurveya, it ta neceeaary that all the accuracy of the field data be 
utilized. Aa will be aeen from the data prevented here, even the full three or four place accuracy of field data 
often allows a fairly wide range of poaalble solutions. The obvioua disadvantage of direct Interpretation la the 
large amount of numerical work Involved as compared with matching apparent resistivity curves. 


The physical situation, aa mentioned above, consists of the halt-space l » o (circular-cylindrical 
coordinates) made up of a given finite number of strata parallel to the surface i * o . Each stratum has thickness 
d L and is homogeneous und isotropic with resistivity (1 • 1. 2, 3 --- n, d n - ) (Fig. 1). On the surface of 

the half space, at / •<?, *••». la a point source of D.C. current, the sink being the hemispherical shell of radius 

R • (r 2 ♦ 1/2 ,a> xhe quantity analyzed 


So^rc 


d, 


is the potential as a function of distance along the 
surface from the point current source. 



r k U 

1 Fi*u »•# 1 

cc 

In any given field situation, the validity of the assumptions of parallel horizontal layering, homogeniety and 
isotropy,if not obvious, can be tested. There are certainly areas where these hold. The assumption of some 
given number of layers is not much of a limitation, since an excessive number can be assumed. The analysis 
should then give to the extra layers either resistivities equal to those of adjacent layers, or zero thicknesses. 
The determination of the effective (f (r) from apparent resistivity values is a simple matter, and is described 
by Mooney and Wetzel <1955), 


Mathematical Analysis 

In the case, as here, that resistivity is a function only of depth . the potential on the surface of the 
half-space can be expressed as (Stefanesco, 1930) 

111 <p(>) • J Irfl) $(M Ul 

where k( X) is a function determined by the way in which the resistivity varies with * , and X , for the pur¬ 
poses of this description, can be called a parameter of integration. The function K X), is called the Slichter 
kernel. For some special kinds of resistivity-depth functions, ^ ( LI can be analytically determined in terms of 
the parameters of the function (Stefanesco, 1930; Slichter. 1933; Langcr, 1933, 1936; Pckeris, 1940; Sunde, 1949; 
King, 1933), This is the case for stepwise variations, such as ours, and for real or imaginary exponential varia¬ 
tions. among others. A systemic representation of t„( i ), the kernel for r layers, appears in Sunde (1949, p. 54). 


( 3 ) 


l J / Alt hf( 


where 


Mu- 


s s 


f'.fc 




(►« ~i)tm 


f -. * e 






Thu expression is obviously non-linear in all parameters. 

Fortunately, equation (1) is a Hankel transform, and can be Inverted to give 

<3) it (a) * 1 f f(y) J 0 

Thls reUtion permits the evaluation of a Slichter kernel from a set of field measurements. 
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Numerical Analysis 

The procedure followed, then, is to assume that * layers exiet and that, therefore, the field kernel can be 
fitted by a kernel k n ( A. ). with the values of the parameters and yet to be determined. The field kernel 
k ( A j) is then evaluated for m values of Aj? 

m > 2 n-1 
J ■ 1 . 2, ---, m 

We then have m non linear algebraic equations in (2n l) unknowns, and the remaining problem is purely numerical. 
General methods for handling this problem are described in the literature (see. for example, von Sanden. 1923, or 
Householder, 1953). The method first used, and described here, is a combined Newton-Least Square solution. 


Let the field kernel k ( X ) calculated for ^ be k ( A,) 
and with initial estimates. X f j\ of the parameters sought. . 

< ( '?)-*. - C' J - K t - 1§ J 


. Assuming an n-layered structure, 


lf„ (\ Jfy ) - K", ' - K 


O 


Here the X*re the P and d., with 1, 2,-, 2 n-1. Likewise, for a second approximation to the parameters, 

•v. < - ,, - <? 

a*?- C-1” - 

The further assumption is made that the kernel can be expanded in a Taylor Series about the value from any 
estimate, for example the first, l.e. 






In particular 


K . k M + l'[ &i) (vO * Cx r“’) 1 

/ 1 A A 


4 cross terms and higher 
derivatives 


♦ cross terms and higher ] 
derivatives 


The Newton Method of first approximation then states 


l l /»> ■ X g 


M 


It is seen that this is an approximation which becomes more nearly exact as the second and higher derivatives, 
and cross derivatives, diminish in magnitude relative to the first derivatives, and that the statement Is exact only 
tn the case that the function is linear in its variables. Such is definitely not the case with the Slichter kernel, as 
is seen from equation 2. However, within a very small range even exponentials can be approximated by straight 
lines, and thus the implicit assumption is that the first estimates X l J' ot the parameters are reasonably close 
to the true values Xj • 


•See the Appendix . or Mooney and Wetzel (1955). The latter is especially well suited to hand calculations. 
••Actually, the values of f and f arc known from the asymptotic values of the apparent resistivity curve, 
and there are only <2n-3) unknowns. 
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Then 


A * (,J 

A * i 


A i in , l 1 , a < 

Ak \ t -. \ ' **. v;' 

The standard Newton Method is to evaluate ^ for <2 nl) values of A. . and to calculate the x3 l, by setting 


1 1 2 t 

l-l 


i( , « *; a 




A t , 


. Co; 


and one has (2 n l) linear equations Ul <2n 1 ) unknowns, the Ay . which can be solved exactly, giving 

■X/ = X *' ♦ A Xf 

[f no errors were present in the field data, and none were introduced in the (numerical! integration of K,. 
then it i ould be expected that kn ( T- ) would exactly equal K ; . However, even barring departures of the 1 

geology (rom that assumed, errors do exist in all field data. Readings of voltage and current cannot be exact 
because of the limited accuracy to which meters can be read. Add to this the usual vagaries of field equipment, 
and it < an probably be said that the K, integrated from the best field data contain random errors in the third 
significant figure, and occasionally even in the second significant figure. As was discussed in the introduction, 
and as documented by the results, at times even the fifth or sixth place can determine changes in a solution that 
may be important when planning a drilling program. Thus, it is desirable to be able to eliminate the contribution 
of random errors, and this is the purpose of the Least Square modification. This formulation minimizes the function 

£ ( _ ( iS ‘ i 

t a. y 

l = I 

with respect to the A Xg . vU: 


4 1 

L l A c A ”) £- ) i A 1 ' 11 * ) = O 


for all P = l 2, ... 2 n-l 


but. from above 


*T~ 7 > ■ V * 

0*1 * g 


A X (St ) 

A x. 


so that 
2* -1 


r 

f 1 


i t- 




; Y <’> 
* / 


for each value of f . This is. again, a set of (2 n-l) linear 
in matrix notation, these equations are 


f /'OV' 

1r, ■ Tx,'x“> 


equations in (2 n-l) unknowns, the A Xy. Summarized 


<4) /o 


a " 1 r 

*e, * 2- 

J a *1 


Vffa ) 

2*/' 


■LfO 


evaluated at X X 


L f ' 1 v" 

b - 2Z * 


x is ^ JX / X irJ 
V * , 


. 1 *) 




The question arises as to the behavior of the solution at minimum error, when this minimum value is 
greater than zero, and when, in parttculai, Uw largest portion of the error is due to a single point. Reference to 
equation 4 shows what will occur in this situation. The partial derivatives change very little over a small range 
of variation of the parameters, whereas the errors. . change relatively large amounts. When the kernel is too 
far from the erratic point, the large from that one point will dictate both the magnitudes and directions of 

parameter changes. As the error due to the erratic point is decreased, the errors due to ine other points (which 
dictate parameter changes in the opposite direction) are increased until the net predicted changes are nil. 

Computation Procedure 

The foregoing was a rather complete statement of the problem and the attempted solution. A program was 
written for the Whirlwind 1 computer to handle the three-layer case. As mentioned in a footnote, in the practical 
case both ^ and £ arc known, hence the unknowns remaining are d^, and d^. The input data were: 

a) the number and values of X 

b) the values P and P , and the estimates of , d, and d„ 

h ij U 1 2 

and c) the 

For most cases, ten values of were used, in the range from A ■ 1.0 to A‘ 0 . 006 . All resistivities were 
reduced to • 1.0, and ^ 2 . dj and d 2 were estimated from the kernel behavior or from apparent resistivities. 

Theoretical three and four-layer kernels and field kernels were analyzed. Several hundred theoretical 
kernels were calculated, with seven to eight place accuracy, on Whirlwind. Field kernels were integrated on 
Whirlwind, as described in the Appendix. 

The machine calculated all the elements of the matrix, solved the matrix using Crout's method (Hilde¬ 
brand. 1952, p. 503), determined the new values of the parameters, and returned to the beginning. It continued 
iterating until the error 22 & 1 began to make random oscillations, indicating no more possible improvement, or 
until the solution blew up. indicating one of several situations, to be discussed in the following sections. 

Several practical computational difficulties should be mentioned. First, the kernel for the three layer 
case is (using the notation of Sunde, 1949, p. 55) 


V s u) - (1 * Ax* 



- 0? * r. J 


« ( J - /'zj ft) /( 1 - /'ll ) 

fa 

* (s-fi j 



and 

t, * c 


T r 

v ■* e 

The algebra 

was simplified by using throughout 

f } n; 

1. l/Cl + l'■,(>)) - 'h * fa T 

instead of 

k,Cx) . 


Second, as will be seen from the results, the parameter changes calculated, especially at the beginning of the 
calculations, are almost always far too large. This necessitates an empirically determined parameter change 
multiplier, M. to prevent the solution from oscillating wildly. The reason foi the overshoot, of course, lies in 
the departure of the kernel from the linear behavior postulated for it (see Hildebrand, 1949, p. 364). Second 
derivatives calculated were, indeed, quite large. Th* valut-s settled upon fui V» were ^ beginning of » 
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calculation. Increasing to 1.0 near the end. if a solution begun to diverge after an increase of M. the neat smaller 
value of M was brought back. 

The Crout solution was used because it is economical of machine storage. It docs not take advantage of 
the symmetry of the matria, and thus would not be used in the case of more variables. Also, the accuracy is 
especially sensitive to the range of variation of the diagonal elements, and it was found quite important *n some 
cases that each row be normalized by this element. 

Another suggestion, ulthough it was not followed in these calculations, is that once the solutions atop changing 
rapidly the derivatives need not be calculated at euch iteration, but rather, the same values may be used for sev 
eral Iterations. 

It is Important to note that the kernels for \ery large and very small values of ) need not be considered, 
since these values wtll be monopolized by -'’and f* , respectively, which are treated as knowns. 

*3 


Results 


Analysis was attempted of 22 three layer cases, 3 four layer cases and 5 field cases. 

The three layer kernels were all known to seven place accuracy. However, the analyses were done on 
data rounded off to three place accuracy, as would be the case for field data. In addition, several of the same 
cases were tested with more places of accuracy to observe the improvement in the solution. Rounding off the data 
in this fashion is equivalent to Introducing a random error of range * i, ln the last figure. Many of the cases 
chosen for analysis were moderately to very difficult in that either dj was much less than dj, or the resistivity 
contrasts were low. or both. One quite interesting result, which is intuitively appealing. Is that in the case of f* 
much greater or much less than /J there is a functional relationship between f *2 an<1 d 2* When f *2 >> 
the product f j ^2* hut not f g or d 2* can h e found accurately. When / j only 'he ratio can be determined 

to any degree of accuracy. That this should be the case can be seen from the asymptotic behavior of the expres¬ 
sion for kf,. 


In several of the cases presented, this procedure was unable to give a solution which was any better than 
the first estimate. This Is probably due to irregularities on the error "surface", along which the iteration is 
proceeding, and could. In some cases, be eliminated by some other kind of procedure. To check this hypothesis, 
a modified "steepest descent" method was tried. Preliminary results indicate that this will help (see p. UO ). 
More work is being done in this direction. 

It will be noted In the results which follow that some kernels have been analyzed more than once. This 
has been done to investigate the reaction of the analysis to 


a) changing the first estimates 

b) changing the accuracy of the kernel 

and c) changing the size of M, the parameter change multiplier 

With regard to a), it wus found that there exists a certain range, about true value, for each of the variables, within 
which the solution will be substantially the same regardless of the first estimates. Qualitative investigation of 
the size of this range leads ue first to the question of uniqueness of the solution. It is the author' a belief, although 
it has not beenjroven, that for completely accurate data only one physically allowable solution exists (excluding 
degeneracy, of course). That this is so for the three-layered case seems intuitively apparent, but it is not so 
obvious for the many-layered rases. Convergence of lids analysis from any starting point whatever, as stated 
above, would require a smooth error "surface" (the multi-dimensional surface generated by the curves traced 
out by the .rror when varying one parameter at a time). The amoother (t.e., eloaer to ftr.t order) the aurf.ee, 
the larger the allowable M and the more rapid the convergence. Inflections or aecontiary minima tn the surface, 
make machine solution of the problem much more difficult. 

Increasing the accuracy of the kernel data ln everv r»«. Improved th. quality o( the aolutton, as expected. 
The amount of improvement fir hi n ui caeca 8. 13, 14. 17 and 19. 
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At any step in the solution there is some optimum value of M which will give the largest decrease in 
error. To find this value, would require calculating the errors obtained using several different values of M and 
interpolating. Since the error calculation la the most time-consuming portion of each iteration, it is seen that 
•uch a procedure would be impractical. Instead, it was found empirically that a small M was required at the 
beginning of the solution, but that after about 15 iterations the process had stabilized sufficiently to allow much 
larger values of M to be used. The normal usage was M • .05 for the first 15 steps, ,5 for the next 15, and then 
M l. ln the event that a step gave a large increase of error, provision was made to bring back the result of 
the previous step and the next smaller value of M. 

As a test of the stability of the procedure, three four-layer kernels were analyzed with the three-layer 
analysis. U is seen that the solutions arc reasonable approximations to the true situations and that the final 
errors do not get as small as with the three-layer kernels. Also, converge is slower than Is the esse with three- 
layer kernels. The same general tendencies might be expected of field kernel analyses. A rather unexpected 
result is that this analysis (in these three cases) tends to ignore the third of the four layers. 

The ultimate test of the method, of course, is whal it will do with field data. Five such cases were tested. 
Figures 2, 3. and 4 show the apparent resistivities, integrated potentials, and kernels. 

Case A did not satisfy the assumptions in two respects. First, the > fC+J&id not reach a constant value, 
and had to be smoothed artificially, and second, lateral variation existed as shown by the apparent resistivity 
. urve (see Fig. 2) and by another taken in a direction perpendicular to that of the one Bhown. Drilling showed a 
good conductor at 6.2 units of depth and below. Note that the solution is degenerate. 

Case B appeared to be a nearly ideal three-layer case. Drilling showed a good conductor at 5.0 units 
depth and below. 

Case C from the apparent resistivity curves, would seem to be at least a four-layer case. No drill holes 
existed in the area. 

Case D from the apparent resistivity, is very nearly a three layer case. The fit obtained is exceptional, 
but no drill holes exist. 

Case E appears to be at least a four-layer case, however, the solution obtained degenerates to a two-layered 
case. No drill holes exist. 


Case 


A 

B 

C 

D 




r. 



THREE 

-LAYER 

ANALYSIS 

OF FIELD 

KERNELS 



First Estimate 

Stillltlol’ 

~7T 

75 

f 3 

£ s 

A 

ft 

ft 

j* e< * 

d i 

d 2 


A ‘ 

1-1 

2l 

ft 


_ 

1 . 

2.0 

1.48 


1 . 

1.53 1.48 


.15 

7.0 


8.5x10’* 

.145 

271.3 


2.4x10’* 

1. 

.6 

.2 


1. 

.641 

.2 


.2 

4. 


1.7x10’* 

.152 

2.68 


1.2x10-* 

1. 

.4 

.21 


u 

.252 

.21 


.5 

4. 


3.7x10’* 

.271 

5.72 


1,7x10’* 

1. 

.75 

.175 


1. 

.708 

.175 


.1 

2. 


1.2x10* 

.243 

.882 


s.ixio-* 

1. 

.2 

.147 


1. 

1.00 

.147 


1 - 2 

5. 


1.2x10 * 

-1.70 

♦2 29 


2.1x10’* 
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THREE LAYEK ANALYSES OK THREE LAYER KERNELS 




C»w 

18a 

16 b 

17a 

17b 

18a* 

18b 

19a 

19b 

30 

21 

22 


THREE LAYER ANALYSIS OF THREE LAYER KERNELS (Contd) 


True Values 

Solution 

First Estimate 


Ft 

Ft 

F» 

Pi* 2 

Fi 

7t 

A 


F 2^2 

A 

A 

A 


I 

d i 

d a 


VA 

d > 

d 2 



V/a 

d i 

d a 

d 3 



1 

.75 

.1 

.075 

1 

3.3 

.1 


.429 

1 

5 

.1 



1 

.1 



.802 

.130 


2x10-3 


.6 

.1 


3*10*5 

3 

1 

.75 

.1 

.075 

1 

1.87 

.1 


.290 

1 

.5 

.1 



1 

.1 



.651 

.15ft 


10-3 


.8 

.08 


SxlO'* 

3 

1 

.05 

.2 


1 

.171 

.2 



1 

.02 

.2 



1 

.1 


2 

.976 

1.12 


SxlO-® 

6.6 

1.1 

.25 


8xl0‘ 4 


1 

.05 

.1 


.1 

.075 

.2 



1 

.02 

.2 



1 

.1 


2 

.997 

.158 


4x10-® 

2.11 

1.1 

.25 


8xl0‘® 


1 

10 

5 

1 

1 

5.00 

5 

-38. 


8 

5 



1 

.1 



.917 

-7.18 


10-4 


.8 

.5 


4x10-3 

3 

1 

10 

ft 

1 






1 

9 

5 



1 

.1 



no improvement 



1.1 

.11 


5x10'® 

3 

1 

1.5 

1 

.15 

1 

.708 

1 




.9 

1 



1 

.1 



1.39 

.0840 


10‘ 4 


1.2 

.ft 


2x10*5 


1 

1.5 

1 

.15 

1 

1.46 

1 


1.57 

1 

1.2 

.1 



1 

.1 



.996 

.108 


5x10*12 


1.1 

.08 


2x10-5 


1 

.1 

10 


1 

.012 

10 



l 

.05 

10 



1 

.01 


.1 

3.51 

4.08x1 <r* 


8x10-* 

.033 

1.2 

.05 


4x10-2 

3 

1 

100 

1 

1 

1 

2.20 

1 


1.20 

1 

5 

1 



1 

.01 



.904 

.543 


10-5 


.8 

1.6 


5x10-2 


1 

2 

10 


1 

31.4 

10 



1 

3 

10 



1 

.01 


.005 

.976 

• -.003 


SxlO" 5 

.00011 

.9 

.1 


4x10'® 



•Degenerate solution. 


THREE LAYER ANALYSIS OF FOUR-LAYER KERNELS 


True Values 

Solution 


First Estimate 


B* ) 

< 




nr 

Pi ft 

d 2 

Error 


% 

Fj 

Error 


l 

100 

1 

.01 

1 

49.2 .01 


i 

50 

.01 

3 

90 ] 


l 

4 


.977 

2.12 

4x10-5 

.8 

2 


2x10-3 


1 

.01 

1 

100 

1 

.0252 100 


1 

.05 

100 

3 

50 

1 

1 

4 


.987 

2.65 


.8 

1 


.4 


1 

10 

1 

.1 

1 

10.7 .1 


1 

ft 

.1 

3 

84 

1 1 

4 

1 


1.05 

3.63 

10* s 

2 

5 


.07 

1 



•Thia column indicates number of Iterations required for solution, with M constant «.05. 
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APPENDIX 

Numerical Integration of the Slichter Kernel 
The procedure used to integrate 

k< l)* + 

o 

was as follows. The potential data V <pCr\ were given for a finite number of points, r . These ure usually at 
30. 50. or 100 foot Intervals, depending on the electrode spread used. The assumptionT>f three layers dictates 
that beyond a certain distance. r c ( V ( >> ) the potential fall off shall behave like that on a uniform 

medium of resistivity ^. and should therefore vanish as I . Thus beyond r c , Y <pfrJshould be constant. The 
value of this constant is determined in this fashion: 

enO r.,) 


f- 4"' ?.,) - 

The quantities r^. r j ^ j • an<1 are * ne *sured. and 

-- YJ 

was extrapolated to 4? .by eye. Data from 20 equally spaced values of r were used, including r ■ 0 and 
r - r c . The data were replaced by two sets of fifth order orthogonal polynomials (sec Milne, 1949. p. 265. ff.). 
one fitted to the last 15 points, the other to the first 6 points, where changes were more rapid. These then gave 
two series which represented the Y <f(Y) in the intervals named 

Y \* ! j ) z c ° * C,py " c * * " C S p r s P - 

v- * d B , j, r;♦ p;, ♦ 4 p; ^ 

Rearranging by powers of r 


J "1 K r 
l - s.t. if 


% * IS * V 


P t - K - 

The integral becomes 








Iffi)* xf Yr p,(-)J/xv)Jy *’ 


where, of course, it is assumed that r constant beyond 

r * r c . 

if 


k(x)* iz f v x Xk( ri 

(? » 

\f l Yli 
the integral 

/ y ' £ ^ 
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can be evaluated ana llyic ally. via: 

b 

-- i-[ j au+ T s u*-) + jjU-)- j k 
»-i5 (*>■)/► - -ji (x u J * 

yVj.awvt- = ovj +s 

* fl'ZUH + it(j r un * (xt.) . o , ;j ‘ 

\/J. (MJ* * j* [ <IW <^TXW t(x>? (X>)]f 

Btc. 


3. ACADEMIC PROGRAM 


3.1 Institute Couraea and Seminars 
MACHINE-AIDED ANALYSIS 

Sixty-one students attended the electrical engineering course 6.25 entitled Machine-Aided Analysis which 
was offered first semester. 1954. The course, under the supervision of Dr. W. K. LlnvlU. Associate Professor of 
Electrical Engineering, was given jointly by Profs. Linvlll. R.C. Booton. and C. W. Adams. 

Both analog and digital computational techniques were studied, as well as the numerical analysis and the 
solution of engineering problems by numerical methods. Although the emphasis was on typical engineering 
problems, some business problems were taken up also. 

As practice in analog computational techniques, the facilities of REAC (Reeves Electronic Analog Computer) 
of the Servomechanisms Laboratory at MIT were used by the students for the solution of an ordinary differential 
equation. 


Thr £ were calculated by the machine using, at most, thirty terms of the series. It was found Important 
that they be quite accurate In absolute value, since some of the coefficients f and k. arc very large. If a 
standard set of values of > were decided upon. It would have been advantageous to tabulate the values of the 
functions. 


Then 


In connection with the study of digital computational techniques, a problem (see Problem #221 in Part 11 
of the present report) was programmed by each student for TAC, the summer session TThrec Address Computer. 
The facilities of the Digital Computer Laboratory were used by the students for the preparation of their routines, 
l ime was made available for each student to be present at the running of his routine, for the correction of mis¬ 
takes, and for the re running of the routine. 



INTRODUCTION TO DIGITAL COMPUTER CODING AND LOGIC 

Sixty-five students are enrolled in course 6.535, Introduction to Digital Computer Coding and Logic. 
being given second semester, 1955. by Mr. D. N. Arden. Staff member of the Scientific and Engineering Compu¬ 
tation Group of the Digital Computer Laboratory. The course will survey arithmetical algorithms used by high¬ 
speed digital computers and a discussion will be given concerning representation of numbers using an arbitrary 
base and the conversion of numbers from one base to another. A brief description of the logical design of a 
simplified computer, practice in writing subroutines, and the design of an interpretive code will also be given. 
Illustrations of some Important mathematical techniques In the practical application of high-speed computers will 
also be considered. The students will be given practice in the practical coding of some of these techniques by 
using a simplified computer simulated by the MIT Whirlwind I computer. 

NUMERICAL ANALYSIS 

This year, for tne first time. Professor Hildebrand Is giving a second semester course in Numerical 
Analysis, M412. This semester the content of the course is made up of the following topics: Least-squares 
approximation, smoothing of data, quadrature formulas of Gauss and Chcbyshev types, harmonic analysis, exponen¬ 
tial and trigonometric approximation, determination of periodicities, Chehyshcv approximation, continued fraction 
expansions and rational-function approximation, numerical solution of partial differential equations. 

MACHINE METHODS OF COMPUTATION AND NUMERICAL ANALYSIS 

The biweekly seminar for the project reported in Part I of the present report has continued into the second 
semester of the academic year. Its purpose and place in the educational program of the project was fully outlined 
in the past quarterly report under the section headed Group Activities. 

3.2 Group Activities 

NUMERICAL ANALYSIS DISCUSSION GROUP 

A group made up of (1) numerical analysis representatives of the Committee on Machine Methods of 
Computation and (2) mathematicians of the Digital Computer Laboratory have continued to meet informally as was 
reported last time. Some of the discussions have been a continuation of last semester's topics (see previous 
report) and a lew nave centered around U*r Uu»l compute tkc first progress report of tin.- piesent publication. 

Further new ideas or results that the group might contribute to. dtacuas, or criticise will be published in later 
progress reports under the same title. "A Study of the Basic Problem of Numerical 'Analysis Expressed in the 
Language of Computing Machines*’. 
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ACADEMIC PROGRAM 

COULOMB WAVE FUNCTIONS 

The programming of Coulomb Wave Function*, being carried out jointly by some physics department 
member* of the Committee on Mu. him- Methods of Computation, is discussed in the appropriate progress 
report, section 2.2. Part l. 

NUMERICAL ANALYSIS LABORATORY 



Project Whirlwind 


1. REVIEW AND PROBLEM INDEX 


The supervision of the Numerical Analysis Laboratory, which ib open approximately 8 hours a week in 
conjunction with Professor Hildebrand's « ourse M412. and the grading of the homework problems are being 
done by M. Doug Us Mcllroy, Philip M. Phipps, and Anthony Ralston. 


This report covers the specific period of December 26. 1954 to March 20. 1955. During this time 66 prob 
Urns made use of 326.27 hours of the 461 hours of Whirlwind 1 computer time allocated to the Scientific and 
Engineering Computation (S&EC) Group. These problems rover some 15 different fields of applications. The 
results of 22 of the problems have been or will be included in academic theses. Of theBe. 19 represent doctorate 
theses, two master's, and one Electrical Engineering. Thirty seven of the problems have originated from 
research projects sponsored at MIT by the Office of Naval Research. 

Two tables are provided as an index to the problems for which progress reports have been submitted. 

The first table arranges the problems according to the field of application indicating the source of each 
problem and the percent of the WW1 machine time consumed. The second table attempts to arrange the 
reports according to the principal mathematical problem involved in each. In each table letters have been 
added to the problem number to indicate whether the problem is for academic credit and whether the prob¬ 
lem is sponsored. 

It is interesting to note that no programmer has reported difficulties due to machine malfunction or 
mistakes <n service routines. In particular. 97% of the assigned machine time was usable. 

Even though no major modifications were introduced into the comprehensive system of service routines, 
the development of new coding techniques by the S&EC Group was extended by the development of translation 
programs fur MIT's Numerically Controlled Milling Machine and for the use of members of the Servomechanisms 
Laboratory in coding for the Univac Scientific 1103 computer. 
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PROBLEM INDEX 


Field 

Dear ripUoo 

Problem 

Number 

%WW1 

Oowrca 

Aeronautical 

Low a*pert ratio flutter 

ITT. C 

.14 

MIT 

Engineering 

Blast response of aircraft 

119. D 

4.00 

MIT 


ttatu upter rotor etability 

Ml. C 

.70 

MIT 


Transient reeponae of aircraft atructurea to aerodynamic beating 

334. C 

1.07 

MIT 

Chemical 

Si* . omportent diatUlation variable enthalpy and aquUlbrium data 

•130. B.N 

.34 

MIT 

Engineering 

Transient effects in diatUlation 

•147. D 

.93 

MIT 


Reactor runaeay prevention 

•131. C 

1.50 

MIT 


Traaeianta in diaUIUUan columna 

941. B.N 

1.11 

MIT 

CleU 

Reapunae of a five story frame building under dynamic loading 

•303. C 

.30 

MIT 

Engineering 

Dynamic analyala of bridgee 

330. C 

1.74 

MIT 


Analyala of two story eteni frame building 

•133. N 

.39 

MIT 

CbMtrlral 

Dynamics and control of packed distillation columns 

>51. B 

.55 

MIT 

Engineering 

Geology and 

Geophysical det* analyala 

IU4. C 

2.34 

MIT 

Oeophyalre 

Interpretation of earlli surface resistivity measurements 

•133. C 

3.33 

MIT 


Dispersion curves fo» seismic eaves 

•313. C 

1.47 

MIT 

Instrumentation 

An uiterpreuve program tu accept mathematical symbols 

ioa. c 

.00 

MIT 

LolKir story 

Servo reapunae to a loaine pulse 

III. c 

.44 

MIT 


Ault* orrvlatum function of auhnutted data 

337. C 

.IS 

MIT 


Ouldanre end rontrol 

330. C 

3.05 

MIT 


Data reduction fur X-l fire control 

344. C 

.51 

MIT 


Flight path of on aircraft during pullup 

IS 3. C 

.IS 

MIT 

Hydrodynamics Lab 

Investigation of turbulent floe 

333. B.N 

.34 

MIT 

Lincoln 

Tracking reeponae character 1st lea of the human operator 

IM. B.L 

09 

MIT 

Laboratory 

Eigenvalue problem (or propagation of E.M. eaves 

103. L 

1.37 

MIT 

Mathematics 

Number uf atructurea of relations on find* set 

•143. N 

.!• 

MIT 

□apartment 

Mr. tunic*1 

Flow of cumpreaaible fluids taerotiiarnvapreaaur) 

130. D 

1.70 

MIT 

Engineering 

Laminar boundary layer of n steady, compressible Oou in 
the entrance rvgton of a tube 

100. c 

LOO 

MIT 


Evaluation of difference diffusion equallun 

•330. A 

.93 

MIT 

Mete urology 

hynoptlc Climatology 

153. 0 

6.10 

MIT 


Distribution uf guatmeaa In tha rraa stroonphtre 

ISO. C 

.13 

MIT 


Computations of ttw fields uf vertical velocity and horizontal 
divergence 

134. N 

9.30 

MIT 


Surface pressure prediction 

347. C 

1.05 

MIT 

Physic* 

('oulomb uavr fmution* 

•123. B.N 

3.00 

MIT 

Department 

Self-conelatent molecular orbitals 

144. N 

3.90 

MIT 


determination of phaae shift* from eapenmanul cross - sections 

103. N 

.07 

MIT 


Overlap integrals of molecuUr end crystal physics 

•171. B.N 

9.7| 

MIT 


Servomechanisms 

Laboratory 


School of Indus ti la I 
Management 
Dynamic Analysis 4 
Control Laboratory 
1-ob. for Nuclear 
Icbnr* 

RoHtrrh Lab. for 
Electronics 
Mtsi ellaneinis 


An augmented ylim • save mat bod aa applied to ooduim 

Study of tha ammuiib molecule 

B«< hang* integral* between raaJ Slalar OrbiUia 

Varutton -perturbation of atomic wave function and energies 

Transformation of integrals fur diatomic molecules 

Neutron-dautaran *. stirring 

Atomic utlagrata 

Eigenvalues for a spheroidal square well 

Self •onmatant calrulatlon of noctaar maaa danalty 

Theory of neutron reactions 

Enaro levels of dUtomlr Hyde Ida a 

Evaluation of two-ceuter mulm-ulor integrals 

Application of augmented plana e»r* method to chromium rrystaI 

DaU reduction prugram. polynomial fitting 

tfubrouttns tor tha numerically controlled milling marhuw 

TranaUtton program for tha nutnariceUy <-.enrolled milling 


Flight in tar captor control 

Dynamic anslvsu of an aircraft interceptor 

Enargy lamia In a apbaroldal puumtial 

Plant survaya for control ayatama by atattatlral n 
f'ryaUl filters 

<'oni|«rehrnalve ayatam of aarvlra routtnaa 
Spar Ml problems (staff training, etc.) 

Library of subroutines 
Into at inn I motility 


Coin par Ison of simple* and r« lu*t ion methods in linear 
programming 

Courae 0.35 Machine aided analyala 
WWI EM A 1109 translation program 


9.19 Univ. of Chic ago 
.50 MIT 

.IT MIT 

l.OT NOT 

.M MIT 

l.TI MIT 

9.50 MIT 

.03 MIT 

.!• MIT 

.15 MIT 

.IT MIT 

5.11 MIT 

.11 MIT 


MIT 

MIT 

MIT 

MIT 

Maaa.Mem. 

Hospitals 

MIT 


Mathematical Pr oblem 

1. Matrla algebra and equations 

Mairli multiplication, addition, dlagonaluatlun 
Matru equation 
Eigenvalue e 

Zeroes of a matrix equation 
Minimise an analytic equation 
Urthonormaluatton 
Eigenvalue a 

3. Ordinary differential equation* 

General ayatam 

Seven nonlinear first order 

Set of nonlinear fir at order equations 

Two second order 

System 

five nonlinear second order 
Second order linear 
Second order nonlinear 

Two linear second order etth periodic coefficient* 

Second order linear 

Set of flret order equations 

Set of nonlinear flret order equations 

Wave equation 

Set of thirteen first order equations 

Nonlinear second order 

Set of eight flret order equations 

Five aquations with 30 sets of initial condition* 

3. Partial differential equations 

Sr lu-oedinger* a equatlun 
Diffusion equation 

Schroedlnger' a equation in aph. coord, etth non 
separable potential 
Second order parabolic 
Flret order linear set 

4. Integration 

Integral evaluation 

Overlap integrals 

Calculation of transfer (unction* 

Autocorrelation and Fourier transform 
Autocorrelation and Fourier transform 
Integral evaluation 
Auto and cross correlation 

Stationary point of • variations! 

Integral transformation 
Autocorrelation 
Hartree- Fork equations 
Overlap Integrals 
Overlap integrals 

5. Statistics 

Multiple time aeriee 

Calculation of the coefficient* of a multiple 
regression system 
Mr gtmelon analyala 
Numerical prediction 

4. Tranecendeutal equations 

Curve fitting 
System of 47 equation* 

7. DaU He due Uon 

General daU reduction 
Mine distribution 
*. Group Theory 

Generation of projection operator* 

Noniaomorphlc relations un a finite set 

• 9. Complex algebra 

Cample* root* end function evaluation 

10. Linear Programming 

Linear programming 


PROBLEM INDEX 


Haste operations 

('rout'a method 

Iteration and dlagonaUsatlon 

Interpolation 

Basic arithmetic operations 
Schmutt process 
Dtagunalixatlon 

GUI's modified fourth-enter Rungv-Kwtta 
Fourth -order Bunge Kutla 
Second-order Hunge KutU 
Extrapolation of difference* 

GUI* a method 
Difference equation 
Pnaer aeries 
Hcctongulur integration 
Fourth order Kungr hull* 

Difference equation 

OUT a method 

Second order Hunge KutU 

Mllna predictor -corrector formula 

Gill* a method 

Finite differences 

Gill' a method 

GUI* a method 

Summation of aerie* 

Racurrence formula 

Matching save functions at the boundary 

Explicit finite differs ra sa 
Finite dtflercnr.ea 


Trapexoldal rule 
Evaluation of analytic forma 
Sunpena' a rule 
Simpeon' a rule 
Simpson' a rule 
Gaussian quadratic 
Simpson' a rule 
Simpson's rule 
Algebraic recursion formula 
Simpson’s rule 
Trapeaoidal plus variational 
Evaluation of analytic form* 
Evaluation of analytic forma 

Prediction by linear operator* 
Inner product* 

Modified Doolittle method 
Quadratic approximation 


Iteration using least a 
Step-by step 


Polynomial fitting, etc. 
Arithmetic operations 


Mai bins genera Uon 
Direct evaluation 
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Table 2-1 Current 1’rohlcmn Arranged According to Field of Application 
(• MIT Project on Machine Methods of Computation) 


Table 2-11 Current Problema Arrantfrd According to thr Mathernatim Involved 
(•MIT Project on MachUut Methods of Computation) 















APPROVED FOR PUBLIC RELEASE. CASE 06-1104 


a. WHIRL.WIND CODING AND APPLICATIONS 




WHIRLWIND CODING AND APPLICATIONS 


2.1 Introduction 

Progress reports a* submitted by the various programmers are presented in numerical order In 
Section 2.2. Since thia summary report presents the combined efforts of DIC Project* 6345 and 6915, reports 
on problems undertaken by members of the Machine Methods of Computation Group have been omitted from 
Section 2.2 of Part U to avoid duplication of Part I. Suitable cross reference has been included in Section 2.2 of 
Part tl for completeness. 

Of the 66 problems described below, twenty-si* (219. 223, 224. 233-239, 241-245, 247, 249-252. 256. 258, 
260. 262. 263. and 285) represent new problems that arc being described for the first time. Ten problems (123. 
130. 187. 183. 186. 189, 211, 222. 254, and 249) have been completed. 


Tables 2-1 and 2-11 have been set up to provide the reader with a convenient Index to various Interesting 
asperto of the problems. Table 2-1 lists the problems according to their fields of application and indicates the 
source of each problem and the per cent of the SfcEC Group' a WWI problem time (326.27 hours) used by each. 
The mathematical problems and procedures involved in the various current problems are tabulated separately 
In Table 2-U. In both tables, problem numbers prefixed by asterisks represent work being performed by mem¬ 
bers of the Machine Methods of Computation Group. 

Letters have been added to the problem numbers to Indicate whether the problem is for academic credit 
and whether it is sponsored. The letters have the following significance: 


A Implies the problem is NOT for academic credit. Is UNsponsored. 

B implies the problem IS for academic credit, is UNsponsored. 

C implies the problem is NOT for academic credit. IS sponsored. 

D implies the problem JS for academic credit, IS sponsored. 

N implies the problem Is sponsored by the ONR. 

L implies the problem is sponsored by Lincoln Laboratory. 

The absence of a letter indicates that the problem originated within the S4EC Group. 

No major changes have been introduced into the comprehensive system of service routines (Problem 100). 
However, i ertain conveniences such as the inclusion of program annotations on the Flexowriter tape, optional 
suppression of the floating address table, use of automatic scope output requests, and improvements In the 
logging procedures have been incorporated. 


„ ti fT'ile C ,r° UP haV ® b ” en “ ltlvC ^ the devel °P n '««t of special pseudo codes and translation 
™ ," 1 CM MIT Num « r *cally Controlled Milling Machine (Problem 250) and the Univac Scientific Model 1103 

(Problem 256). 


2.2 Problems Being Solved 

100. COMPREHENSIVE SYSTEM OF SERVICE ROUTINES 

The compr.han.lve system of service routine, ha. been developed by the Scientific and Engineering 
Computation (S4EC) Group to atmplify the process of coding for WWI. The system now in use called CS U 
was described in Summary Reports No. 38, 37. and 38. 


tau provides for: 


. ... 77 ; .. r4C *°wr»er-cooed perforated paper tapes; 2 ) the logging of 

each problem on film and paper tape for subsequent processing. 3) assembly during read-in of a suitable set of 
instructions including interpret,ve programmed ar.thmetlc (with an optional floating-point number system), up 
to several hundred cycle counter. (B-boxes). output routine., error detection, and automatic post-mortem.. 

,nd 0 r^T^^Zr° rmany l r ded r h mn,!mon * c °Peratlons, symbolic addresses, relative addresses, progrsr 

^red on .rm .ro;”^ T Ud r " ‘ U ’ d SP * CU1 C ” n,ro1 Wurd *- Thc routines -re permLotly 

stored on magnetic tape and are transferred to s magnetic drum for automatic selection during resdHn. 

mem. made d “ rln ‘ <h “ q “ r “ r ^ * CS a H ™' Improve - 


Comment Word s 

The CS conversion program has been modified so that comments can be typed on a Flexowriter tape along 
with the words in the coded program. A comment must be preceded by s vertical bar and terminated by a 
carnage return. Comments are ignored by the conversion program. For example, 

si 78 | Select magnetic tape unit 1. 

ca al 

rc | Record a word. 

Kcvislon of Automatic Output Routines 

The automatic output routines DIH, DOB, iDIB and IDOB have been rewritten and shortened in the process 
The programs have been modified to include recording on and reading from the buffer drum. 

Flad Table Suppression 

Printing of a floating address (flad) table as a result of a CS conversion can now be suppressed st the 
option of the programmer. This can be done either by using the word cj. In a director tape or by setting up a 
number in the selector panel and using the manual read-in mode. 

Modification of the Logging Program 

The logging program has been modified so that time entries are made in the log both at the beginning of 
a magnetic tape search for a utility program and at the termination of the read-in of the utility program. The 
resulting time increment will be subtracted from the total run-time for a problem by the automatic biweekly 
program. 

Mo dification of the Stop Instruction 

The STOP instruction has been modified to make a time entry in the log whenever one Is executed. The 
exact time at which a program was completed is thus known. 

Modification of the Title Recording Program 

Tape titles have been recorded for delayed printing by the utility control program whenever CS Flexo (fc) 
tapes or post-mortem request (fp) tapes are read in. The program has been modified so that the date and time 
will be recorded along with the title. 

The reference date and time are recorded on the buffer drum by the operators at the beginning of the 
computer operating period. An automatic indexing of the date occurs if the operating period extends beyond 
2400 hours. 

Automatic Biweekly Program 

An automatic biweekly program has been written which will process all logging tapes produced during a 
biweekly period and produce a typed summary of computer operation for that period. 

The typed summary consists of: 

1. time used by each problem 

2. lost time due to machine malfunction 

3. accumulation of all the time used by all the problems 

4. time used for checking auxiliary equipment 

5. number of problems run 

6. number of programs run 

7. unused Ume (time between run*). 
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1 'Hf program will automatically deduct time due to operator 1 a error**, tape room errors, etc. It will 
aUo deduct time used for searching purposes, for instance searching for CS II or the scope post-mortem on 
magnetic tope. Thla program will print on the direct Flexowrlter the date and initial time each logging tape 

started. 

The use of logging tapes has reduced the time needed to prepare the S4KC biweekly report from more 
than 15 h^ur - to about 20 minutes. This new procedure has eliminated most of the day-by-day clerical work 
that was needed to prepare the records pertinent to the biweekly statistics. 

The time used in preparing the biweekly program is logged under problem 100. The logging routines wil] 
continue to be studied with regard to increasing their flexibility and functions. In particular it is planned to 
include in these routines the actual billing for non-ONR, non-Lincoln sponsored projects. 

Automatic Scope Output Requests 


algebraically coded tape, in order to screen out details and leave the way open for a reasonably straight¬ 
forward logical scan. The logical scan is itself only partially complete, however. In addition, much work 
remains to be done on the detailed routines for many of the special features to be provided in the new system. 

It is currently expected that active work on this problem will resume about the middle of the next quarter. 

120 B. N. THERMODYNAMIC AND DYNAMIC EFFECTS OF WATER INJECTION INTO HIGH- 
TEMPERATURE HIGH-VELOCITY GAS STREAMS 

This problem is connected with the development of a potential gas turbine component, called an 
"acrothermopressor", in which a net rise in stagnation pressure of a hot gas stream is brought about by the 
evaporation of liquid water injected into a high-velocity region of the flow. The concepts underlying the opera¬ 
tion of the aerothermopressor are an outgrowth of comparatively recent work in the field of gas dynamics, and 
it* proposed function in the gas turbine cycle is analogous to that of the condenser in a steam power plant. 


I'he automatic output instruction SOA (tSOA). Scope Output Alphanumerical, was made available to 
programmers using CS 11. The special instructions COLUMN 1COLUMN and FRAME 1FHAME were also 
introduced. A detailed description of these instructions may be found in Memorandum DCL-48. 

106 C. MIT SEISMIC PROJECT 

A. discussed in various previous reports. Problem 106 is concerned with the Investigation of the use of 
statistical analysts techniques to seismic record interpretations, and in particular to the separation of "reflec 
ttuns from background Interference on these records. More complete descriptions of the problem and the 
approaches used are contained in the bt-weekly report of June 15. 1053 (M-3244) and In "Detection of 
Reflections on Seismic Records by Unear Operators". (Wadsworth, Robinson, Brv.tn. and Hurley - - GEOPHYSICS 
Vol. 18. No 3 July 1053). 


I he most recent approach to our problem has been a study of the statistical nature of the interference 
mentioned above and tts relation to the physical mechanisms of generation. The hope ts that this study should 
provide a basis for developing and for evaluating existing procedures of noise removal. 


Computationally we are dealing wtth the funettone of generalised harmonic analysis correlation, spectra 
coherency, linear operations, on discrete multiple time series. Computational details are given In previous sun. 

Estimates of these quantities on actual data are compared with theoretical and experimental 
estimates on truly staUonary series, and interpreted in the light of various types of decomposition! 

arttv" 1 rave 1Z°''' nr '"' d *ra* ” t * mfc ‘ nt ' p, ' r ' n ' " approximate. stationary series, but that the "station- 
nrarlv Tteo t .L h k T 2 XT’ * ,r0m ab "° ,u *'' “■‘‘“‘V <* *h« integrated power spectrum to 

.- ~ 

obtainecTthrough Hte^pecUl^.ollTctton ^ * 

sen, the companies. Occasionally pC»«^XhY^' “ “* 

The programming has been done by S. Suupson D Gnne S Treit.i ,.„a i v i . 

assoc iated with the MIT Department „r Geology and Geophy.ic^’ * ° ( 

108 C. AN INTERPRETIVE PROGRAM 

routine a permitting the tran^tto^^W^To^a'probtom^trt W(Urlwlnd 1; ,e " * " ct of 

sre at present J. H. Utning. Jr. and C. Block of the MIT InatntmenUtton programmer 

.he pj^era resrs? *** ^ . .«*■ **»>« 

, gran. nave even written for me execution of two preliminary scans of an 
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The device consists of a converging nozzle which accelerates the exhaust gases from the turbine into a 
circular duct of varying diameter terminated by a conventional conical diffuser, which recovers the kinetic 
energy of the flow before discharging it to the atmosphere. At the entrance of the duct, special injectors 
deliver minute jets of water which are in turn atomized by the rapidly moving gas stream. 

The changes in state within the aerothermopressor are brought about by the simultaneous thermodynamic 
and dynamic effects of (a) evaporation of the liquid water, (b) momentum and energy interactiona between the 
phaaea. (c) friction, and (d) variations in cross-sectional area of the duct. Under proper circumstances, these 
effects bring about a net rise in stagnation pressure across the device. Further descriptions of this device may 
be found in earlier reports, beginning with Summary Report No. 32. Fourth Quarter, 1052. 

The role of Whirlwind I in the successful development of the aerothermopressor is intimately connected 
with the determination of performance characteristics of the device under all conditions of operation by means 
of a comprehensive one-dimensional analysis of the process. This analysis involves the simultaneous solution 
of seven non-linear, first-order differential equations. 

During the past quarter, a completely revised Whirlwind program treating the aerothermopressor nnaly 
.sis was completed and checked by systematic testing. This program, the writing of which began December 1, 1054 
and required approximately ten weeks, involves the use of 5370 registers. The speed of computation, compared 
with the original program, written in early 1953, has been increased from four to six-fold, the latter figure for 
the situations in which a special wet-bulb temperature computation is required. This has been accomplished 
by elimination of interpretive programming wherever possible. 

This program is directed by a single parameter tape, which will eventually be prepared by filling out 
a standard form, snd therefore represents an automatic computation facility for the aerothermopressor. 

Some of the features of the program are: 

1) Inclusion of five different numerical methods for solving the differential equations. These are (a) 
Euler method, (b) backward differences, (c) second-order Runge-Kutta. (d) fourth-order Runge-Kutta. 
and (e) forward and successive differences. The computations may be directed to change from one 
method to another arbitrarily. 

2) Provision for automatic change of increment (according to an arithmetic progression) at each step 
if desired. 

3) Automatic detection and elimination of oscillations In liquid velocity. 

4) Determination of singular solutions of the equations by a completely automatic iterative procedure. 
(This procedure, which requires about 5 minutes per solution, originally required about 30 minute h 
with the need of intermediate hand computations). Early stages of the iteration use the Euler 
method with an automatic change to fourth-order Runge-Kutta when more accuracy is required. 

It is planned to resume computations on u broad scale. Considerable experimental datu from the 
medium-scale aerothermopressor test facility la now available for corroboration of the analysis. Further, 
several procedures for calculating optimum performance of the aerothermopressor and for studying the effect 

of vanuluios in droplet apfetrum nw*lt exploitation. 
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The ae rot he r mopre s sor development program is being carried out at M.I.T. under the sponsorship of 
•h«* Office of Nava! Research and is directeJ by Pi ofcaaur Ascher H. Shapiro of the Department of Mechanical 
Engineering. The theoretical aspects of the problem treated by Whirlwind I are being carried out by 
Dr. Bruce D. GavrlL 

122 N. COULOMB WAVE FUNCTIONS 

A report on this problem is given in section 2.2 of Part 1. 

ILM N. EARTH RESISTIVITY INTERPRETATION 

A report on this problem is given in section 2.2 of Part I. 

126 C. DATA REDUCTION 

Problem 126 18 a very large data-reduction program for use ut the Servomechanisms Laboratory The 
overall problem is composed of many component sections which have been developed aeparately and are now 
he mg , ombtned into c omplete prototype programs. Descriptions of the various component sections have 
appeared m past quarterly reports. After the development and testing of the prototype Whirlwmd programs 
..re completed, the programs will be re-coded for other, commercially available, large scale computers 
(probably the ERA 110.1. IBM 701, and IBM 704 computers), for use by Interested agencies for actual data 
"•duetto" at other locations. The programs are currently being developed by Douglas T. Ross and Walter Whelan 
th h " of Miss Dorothy A. Hamilton. Servomechanisms Laboratory staff members. This work is 

sponsored by the Air Force Armament Laboratory through DIC Project 7138. 

the „rnIam na .!'!. rC , 0f PraW *T re '” i ‘ rB8 not anl y extreme automatlcity and efficiency in the actual running of 
RJam. but also requires the presence of human operators in the computation loop for the purpose of 
■h ision making and program modification. For this reason extensive use is made of output oscilloscopes for 

w h the ° l T r x 1 un ". nu "“ and manUi “ >"«RRventian registers so that the operator ran cormmmirate 

i . ompu cr. The intent is to allow the human operator to communicate with the computer in terms of 

detTld st“; 6 C °7 U,er *’ rUnmne> and have ,he program leansia^ these ideas inU,T 

, . * . ' !c " ' ***** ° r P ru * ram modification to conform to the human operator's decision The program 

h does this translation and modification is called the Manual Intervention Program (M1V) The . nrrent 
vers,on of the prototype data reduction program is called the Basic Evaluation PrfgTm 

^ium. (in March ninth the agenda included » descritfion^nd d" agencies attended the sympo- 

•" . . .. -i—b*< ... 

is being expanded to include /ne w7ldmunil oTequattona"a^toe^W^^ Pr ° grame - TbB ba8ic Program 

The new logic will be limited only by the drum capacity of the compute"/ 

Arrangments were made with members of the S anH pr ~ » 

utility programs for use on Problem 126 Three routine ^/ UP K« 0r hh0rt tCrm develo P ment of several 

routine, a logging and editing rouUne o'gw record, of .D° m ,i V'y ^ ”’ mplB,Bd - ^>ude * »eope 

■f the Director Tape program. These rou^s m L dUr ‘ ng * “ d “ ™ dlf ‘ Bd *•"*« 

in later quarterly reporte. m thc new P rot <rtype program and will be described 

initiated by thc sponsors of Problem^26 Tdic 7 ^£$) * AT 1103 com P uter - Problem 256 was 

elsewhere in this quarterly report * escri Ption of the problem and current progress appears 
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130 B. N. SIX COMPONENT DISTILLATIONS, VARIABLE ENTHALPY AND EQUILIBRIUM DATA 
SIMULTANEOUS NON LINEAR EQUATIONS 

A report on this problem is given in section 2.2 of Part I. 

131. SPECIAL PROBLEMS (STAFF TRAINLNG, DEMONSTRATIONS, etc) 

This problem has been set aside to account for the WW1 time expended in demonstrating the computer 
to visitors at the Digital Computer Laboratory, in developing routines to be used in these demonstrations, and 
in testing routines written as part of the training of new staff memtxTS. 

During the past 3 months, 13 groups visited the Laboratory. The affiliations of some of the larger groups 
.ire given in Appendix 2. 

A calendar routine was prepared and tested by R. J. Hamlin, MIT Staff member in the School of Industrial 
Management and E. Half fa. of the SAEC Group for use in demonstrations. Given any date from March 1,0001 to 
December 31. 9999. the program calc ulates and prints on the direct typewriter the day of the week on which it 
falls. If the given date should fall on a holiday, the name of the holiday will also be printed. In addition, the 
program calculates the date of Easter Sunday for any given year within the above range. 

132 D. SUBROUTINES FOR THE NUMERICALLY CONTROLLED MILLING MAC JUNE 

At present, milling-machine tape-preparation programs for two different types of cams are being 
developed. One program has apparently been tested successfully and should produce a useful tape during the 
next run. The program for the other cam is nearing completion. Writing and testing of additional subroutines 
is continuing. 

John Runyon of the MIT Servomechanisms Laboratory is using his study of methods of miUrng - machine 
data preparation which Included the development of the subroutine library to facilitate the tape preparation 
process as research for an Electrical Engineer degree. The thesis is expected to be completed this term. 

141. S4EC SUBROUTINE STUDY 

Gill's modification of the classical Runge-Kutta fourth-order method of solution of systems of first- 
order differential equations, described in detail by S. Gill, Proceedings of the Cambridge Philosophical Society, 

47 (1950) 96-108, has been coded and put in subroutine form. The subroutine finds the solution at the end of one 
interval of the system 

5^ ■ f i <*• *V .^n' * 1 - «<1> n 


y t <* 0 * * lv< ' n 


x x Sx„ 


♦ h 


An auxiliary subroutine to evaluate the functions f^ must, of course, be written for each different system. 

A subroutine to find a relative minimum of a function f of n variables has been written. The method is 
to start from an Initial estimate x and proceed in the direction of the negative of the gradient for a distance ^ 
The value X is also estimated to begin with, and on succeeding iterations Ih taken as equal to, half, or twice 
the previous value of A according as the cosine of the angle between the gradients of corresponding values of 
t is between 0 and .9. negative, or between .9 and 1. The process « ontinucs until successive values of f differ 
at most by some preassigned £ . 

A method of economization of power series, by means of Chcbyshcv polynomials has been coded. The 
procedure, descrll^ed more fully (e.g.. by C. Lanczos, Tables of Chebyshev Polynomials, NBS Applied 
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Mathematic* Series. no. 0) is briefly as . a power series « n ia truncated and then expressed in 

term* of Chebyshev polynomial*: a x 


_a 

X X 

n -0 ft . «0 


» t ft T 

n nA k 


X! ( XI ‘nk) T * 

kftO n- k 


where the inefficient* t ^ nifty be founo in the MIS publication mentioned above. Now the T. , may be omitted 
one by one until the uniftrm error so m i ie reaches a desired margin. The polynomial approxfniation may then 
be rearranged In powers of x again giving in general a short polynomial which approximates the original function 
more closely than the power series if cut off after the same number of terms. 

144 N. SELF-CONSISTENT MOLECULAR ORBITAL 

The program represents a mechanization of the Roothaan scheme for solving the many electron self-con- 
*‘Stent field problems within the framework of linear combinations of specified functions. The block diagram and 
details of thr program have been given In a previous Summary Report. There Is no mathematical guarantee that 
the procedure should converge and the current difficulty has seemingly arisen from using input data not possess ir.p 
the proper symmetry properties. It is now felt that the program Itself la operating and further work on it has 
been turned over by Dr. A. J. Meckler to members of the Solid State and Molecular Theory group. 

1SS N. SYNOPTIC CLIMATOLOGY 

The work of the Synoptic Climatology Group of the MIT Meteorology Department over the past quarter has 
been confined mostly to evaluating cloud seeding experiments. The statistical approach used for the past years 
by this project was adopted as the method of analysis. It was of interest to determine how much additional tnfor 
mation was contained in the circulation pattern regarding precipitation over and above the amount usually 
accounted for by control and target relationships. The analysis as far as machine computations are concerned 
has been completed. Evaluation of the results, however, has not been completed. 

162 N. DETERMINATION OF PHASE SHIFTS FROM EXPERIMENTAL CROSS-SECTIONS 

A phase shift analysis is being made by Dr. F. J. Eppllng of the Laboratory of Nuclear Science of the 
elastic scattering of protons by 0 1 over a range of bombarding energies from about .5 Mev to 4.6 Mev. Cross 
sections were measured at eight scattering angles from 168.0° to 90.4°. The cross section i a can be 
expressed a. a function of the scattering angle. 0 and S ♦ , the non-coulomb phase shift of the partial wave 
of orbital angular l and total angular momentum J ■ J I 1/2. 

A program has been written by Mias E. Campbell which will find the S - ’ a that make the sum of the 
squares of the errors between the observed and computed J*'. a minimum. This program has been designed 
to work for any number of phase shifts from 2 to 7 and forTny energy. As soon as the experimental data tas 
been put rnto usable form, the program will be run to determine the number of phase shifts and their values re- 
qulred In the energy range mentiuued above. 

167 B. N. PRODUCTS OF BATCH DISTILLATIONS WITH HOLDUP 
A report on this problem is given in section 2.2 of Part 1. 

172 B. N. OVERLAP INTEGRALS OF MOLECULAR AND CRYSTAL PHYSICS 
A report on this problem la given In section 2.2 of Part I. 

177 C. LOW ASPECT RATIO FLUTTER 

An over-all outline of the problem is given in Summary Report No. 38, Second Quarter. 1954. 
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This involve* the multiplication of five complex matrices by five real matrices and the addition of thr fivr 
resulting matrices. The final complex matrix Is an eighteen by eighteen matrix which results in a system of 
thirty-alx real simultaneous equation*. Coding for the solution of these equations is now in progress. 

This study is being conducted by John R. MsrtucceUl of the M.l.T. Aero-Elastic and Structures Research 
Laboratory. 

183 D. BLAST RESPONSE OF AIRCRAFT 

Since the writing of the last quarterly report, it was found that the Milne's method of numerical Integra 
tu>n may cause instability in the aolullon if the structural discontinuity at buckling is severe. For some para¬ 
metric values, even an additional iteration cycle could not overcome the difficulty. 

As a consequence, the problem wbb reprogrammed by Y. Shulman of the M.l.T. Aeroelustlc and Struc¬ 
tures Research Laboratory using the Rungo-Kutta method. The two simultaneous second-ordrr equations arc 
written into four simultaneous first-order equations to utilize the subroutine available in the SAEC Group Library. 
A test run wss made for a quasi-steady damping case. The results agree with those obtained by an analytical 
solution to five significant figures. 

The peak responses for a given set of parametric values were computed for 300 cases to examine the 
effects of various parameters involved in the problem. The lethal criteria for two fictitious airplanes were 
also established by choosing four values each for altitude, forward velocity, and positive duration of the blast. 

The machine-computing phase of the problem is now completed. 

This work was carried out by H. Lin of the M.I.T. Aeronautical Engineering Department and will be 
Included in his Sc.D. thesis. A modified version of the study is to be published as a U.S.A.F. technical report 
at a later date. 

186 B,L. TRACKING RESPONSE CHARACTERISTICS OF THE HD MAN OPERATOR 

This problem seeks to determine the response characteristics of the human operator as a component in 
a control system. A more detailed description of the problem may be found in Summary Report No. 38. The 
study ia being carried out by J. 1. Elkind of the Lincoln Laboratory at M.l.T. The results will be included in 
his Sc.D. thesis to be submitted to the Electrical Engineering Department. 

The Fourier transforms of several more correlation functions of human operator inputs and responses 
were obtained on WWI. These transforms complete the work to be performed on WWl. The investigation of 
human operator characteristics continues, however, but the power-density spectra required to determine human 
operator characteristics will be computed with the aid of some special purpose analogue equipment which has 
been constructed for this purpose. 

189 C. DISTRIBUTION OF GUSTINESS IN THE FREE ATMOSPHERE 

The statistics of radar weather echoes are related to the statistics of atmosphere turbulence. (Research 
Report 22 A and B. Weather Radar Project). The precipitation is here being used as tracers of the motion. 

The first probability density of the gustiness is obtained from the transform of the square root of the 
auto-correlation function. Six such transforms were computed on Whirlwind using a program devised by Douglas 
Ross under Problem 171 (see Summary Report No. 38). The computations were entirely satisfactory, and for the 
present no more of these are being contemplated. 

193 L. EIGENVALUE PROBLEM FOR PROPAGATION OF ELECTROMAGNETIC WAVES 

This problem was described in Summary Report No. 39, July - Sept., 1954. It arose at Lincoln Labora¬ 
tory in connection with the problem of calculating the electromagnetic radio frequency field radiated by a 
Hertzian dipole to points well beyond the horizon over a perfectly reflecting earth through the lower atmosphere 
idealized to be an inhomogeneous medium with index of refraction decreasing linearly with height. 

Program tapes have been made, based on both the power series and the asymptotic series for the Bessel 
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functions. These arc being used to give numerical results at three different frequencies. It is expected that 
enough results will soon be computed to give n comparison with practical measurements that have been made. 

194 H.N. AN AUGMENTED PLANE WAVE METHOD AS APPLIED TO SODIUM 

The program which generates the wave vet tor groups of vectors of the reciprocal lattice for the body- 
. entered, fa< e centered and cubic lattices has been written and tested. The program will also generate all 
w >» v <• v« . tors equivalent to the reduced wave vector for a given set of neighbors in reciprocal space. Only those 
wave vc* tors which cannot be carried into one another by operations of the wave vector group are generated. 

The projection operators for the various irreducible representations of the crystal types named above are also 
generated on demand of a simply encoded set of numbers. If projection onto an irreducible manifold is desired 
without regard to a particular basis, the routine will generate that projection operator which leads to the least 
number <>f fun« lions in the linear combination which results from the application of the projection operator to the 
function to he symmetrized. 

A program is being written which will take quantities generated by the radial function routine and incor 
poratc them into a tape which will represent in a sense the subspace in which the Hamiltonian matrix is to be 
diagonalized. This tape is then fed Into the matrix diagonalization routine along with the particular symmetry 
de ired and how much of the Buhspa* e represented on the tape is to be used. The way this is done is to have the 
radial wave function routine make a tape which reads into the drum group whose number is the same as that of 
the hand to which the wave functions belong; each successive block on a particular group represents a successive 
neighbor in reeiprot ml space. To generate the tape to he read into the matrix routine a tape is read into high 
speed storage which contains the hands desired, each preceded by a negative number whose magnitude is the 
urd< r of the neighbor in reciprocal space desired in that band. The input Is then a set of coupled numbers. 

Production has begun on sodium and the results obtained are being studied by Mr. M. M. Snffren of the 
Solid State and Molecular Theory Group. 

A program to do the method of symmetrized augmented plane waves has been coded and is heing tested. 

195 C. INTESTINAL MOTILITY 


Ur. .!. T. Farrar of the Gastroenterological Section of the Evans Memorial Hospital is studying the effect 
"f r, *dintion upon th<- motility of the small intestine in the rabbit. The analysis of the records is being performed 
using autocorrelation and Fourier transform, both performed on W.W.l. 

During this quarter, twelve autocorrelations have been completed and the results plotted photographically 
using scope program. Insufficient results are available to draw even tentative conclusions. 

It is planned to perform Fourier transform on all autocorrelation®. The evaluation of irradiation effect on 
motility wtll be based upon differences in: (1) the mean square value. (x 2, , and (2) relative power contribution 
of various frequencies as derived from the frequency spectrum. 

Iflfi. SINGLE ADDRESS COMPUTER 


Sim ,- the 195.1 Summer Session (SS) and SAC programs were written, several changes have been made in 
the locations of the mp-flops in test storage. In order to avoid the restoration of these flip-flops to their original 
addresses when SS or SAC problems are run. SS and SAC are being modified to conform to the new addresses of 
thi flip-flops. The corrected SS has been tested and SAC is now being tested. 


199 N. LAMINAR BOUNDARY 
REGION OF A TUBE 


LAYER OF A STEADY, COMPRESSIBLE FLOW IN THE ENTRANCE 


theorelt a mvel ,1 0,1 ''-‘•'transfer coefficients for supersonic flow of sir In s tube, a 

' 'rncd cM,' bv PrlT T Y r rha " a ‘ 4 " UUl “ ° f ,h « h ' n " n: ‘ r tottery layer in the entrance region has been 
ronUnuilv moment,'.n - 1 of the Mechanical Engineering Department. Partial differential equations of 

mu, T ricToT .c " „y Cr ®' * ere d ' w, °P e<t ,or ,h * ““"‘"or boundary layer. These were then trnnsformec 

condition, at tL tTZll ** ' , ° lV ' ,d f ° r ,peciflc M ‘“ h «-rn»l 


Gill's method is used in the numerical integration of these differential equations. 

Solution of the third set of the differential equations for the case of * onstant viscosity and thermal con 
d:.v tivny has been obtained for adiabatic flow at an entrance Much number of 2.8. However, the convergence of 
the iteration scheme uBed to find the correct initial conditions is very slow ui this case. Work is being done to 
improve this situation. 

The algebraic progrum developed by Dr. J. IL Lulling in problem 108 iu*s» been used in the numerical 
integration of the three sets of differential equations. Recently, a new program using US was prepared and run 
successfully for the solution of the first set of differential equations. The subroutine for the Gill's method was 
also developed by Dr. Lining. This new program is to be used for the solution of the < ase of temperature- 
dependent viscosity and thermal conductivity. 

Solutions of the boundary-layer equations for the case of t onstant viscosity and thermal conductivity are 
now being prepared for publication in some scientific journal, and as a project report. 

201 N. STUDY OF THE AMMONIA MOLECULE 

Due to inconsistencies in the original self-consistent field solution, the self < onalatent field part of the 
ammonia molecule calculation has been redone. In addition, the matrix diagonu Illations mentioned in the last 
progress report have been carried out using the secular equation ruutine, developed by F. J. Corbato under 
problem l72,and the ground state energies of the Ammonia molecule have been calculated. The analysis of the 
solutions is now being carried out by Prof. Harvey Kaplan at the University of Buffalo. This work was • arried 
along by members of the M.I.T. Solid State and Molecular Theory Group. 

203 N. RESPONSE OF A FIVE STORY FRAME BUILDING UNDER DYNAMIC LOADING 

A report on this problem is given in section 2.2.of part 1. 

204 N. EXCHANGE INTEGRALS BETWEEN REAL SLATER ORBITALS 

This problem is being studied by members of the Laboratory of Molecular Spectra and Structure of the 
University of Chicago in cooperation with the MIT Solid State and Molecular Theory Group. A routine was 
developed by P. Merryman of the Chicago group for evaluating two center exchange integrals. A more detailed 
discussion of this routine was given in Summary Report No. 39. The routine is still being tested. 

211 C. SERVO RESPONSE TO A COSINE PULSE 

A description of the original problem may be found in Summary Report No. 39 for tiie 1 hird Quarter. 
1954. In Summary Report No. 40 there is a terminating report for the original problem under thus number as 
well as a description of a related problem. During the past quarter, some work was done on the related prob 
lem, described as follows. 

It is desired to find numerically the maximum response s(t) for some parameter pairs < C , l p / i n >. 

0.01 ± 0.1£ T ll i 10 in the equation 

p n 

(1/ u> 2 ) X -v (2 £ I CO ) X + x ■ G(t), where 
n n 

A/d. 2 . 

2 t (T IT ) sin (2wt/T >1 , 
n p p J 


and G(t) » 0 for t > T p . 


T * 2 T /*»» , x ■ dx/dt, H ■ 
n n 

Git) - 1/2 fl - ros (2 ’Tt'Tp* + 
for OitiTp. 
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Coding of the problem using CS 11 was accomplished by Dr. J. M. Stark of the M.I.T. Instrumentation 
Laboratory. For application of the results obtained, reference is made to pp. 7 04. 705 of INSTRUMENT 
ENGINEERING, Volume II, by Draper. McKay, and Lees, and to Volume III (yet to appear) of the sume work. 


217 N. VARIATION PERTURBATION OF ATOMIC WAVE FUNCTION AND ENERGIES 
A report on this problem Is given in section 2.2 of Part I. 


212 H. N. DISPERSION CURVES FOR SEISMIC WAVES: MULTILAYERED MEDIA 


18 N. TRANSFORMATION OF INTEGRALS FOR DIATOMIC MOLECULES 


The program for calculating dispersion curves of surface seismic waves on multi-luyered media is 
up*'rating tatlsfartorily. As mentioned in Summary Report No. 40, the program was written to handle a five- 
layered medium. However, by resetting one constant, as many layers as desired can be handled. On the aver 
lgc. ime data point for a five layer case requires approximately 1.5 minutes. Reducing or increasing the 
number of layers gives a corresponding reduction or increase in computing time of, roughly, i/4 seconds per 
layer. 

This method is merely a solution for the roots of an algebraic equation, the lowest of which correspond 
to Rayleigh mode propagation. It has been used with no modification to Bolve for the large roots corresponding 
to normal modes of propagation. The method requires twi initial estimates of the phase velocity at each value 
of wave number for which a phase velocity is sought and the final solution is the root nearest to the first esti 
mates. Thug .» Rayleigh or a normal mode velocity will be obtained depending on which is closer to the estimab 

The program has been chec ked with dispersion curves in literature and is being used to calculate new 
results. It has hern derided not to calculate group velocity curves since machine time would be better expended 

on phase velocities. 

An interesting, hut time-consuming, possible application was brought up in a discussion of the problem 
with the head of the Lumont Geological Observatory at Columbia University. It is believed that many unexplained 
t' ults of seismic measurements at sea are due to smooth changes in velocity and density with depth. This 
ituutton might be approximated by a large number of layers with small changes between each. 

This work is being carried out by K. Vozoff of the MIT Geology and Geophysics Group as part of the 

requirements for a Ph.D. thesis. 

215 B. DYNAMIC BEHAVIOR OF INDUSTRIAL PROCESSES 

The ipph> ation of automatic control to industrial processes would be greatly facilitated if the dynamic 
behavior of the process could be determined from measurements made during normal operation of the process, 
without shutting the process down or introducing inputs of a transient or periodic nature. An experimental 
measurement of the impulse response of a linear heat transfer system subject only to random disturbances of 
its operating point has been carried out in the Process Control Laboratory with the cooperation of Professors 
D. P. ( umpbell and L. A. Gould. Input auto-correlation functions and input-output cross-correlation functions 
were computed using programs developed under problem 107 by D. T. Ross. The frequency response computed 
from the results of the correlation study agrees well with that measured on the sume system by conventional 
techniques. 

\ method which avoids the necessity of solving an integral equation to find the impulse response from the 
or relation functions is described in a forthcoming R.L.E. Progress Report, covering the period December 1st 
to Feb. 28. 1955. An autocorrelation und Fourier transform required by this method have been carried out 
using Ross’ a programs. butl.B.M. equipment will be used for the less massive computations which remain to 
be done. 

The results of this experimental work have shown that general purpose digital machines, such as 
Whirlwind I. as well as business ma.-hines of the I.B.M. type. are much more suitable for computing the corre¬ 
lation functions required in the study of industrial processes than are special-purpose analog correlators. The 
digital machines offer higher accuracy than the analog machines, and are free of limitations of available delay- 
time und frequency range. 

Tins work was carried out by S. Margolis of the Statistical Communication Theory Group of R.L.E. under 
the supervision of Professor Y. W. Lee. 


A discussion of this problem has been included by Dr. Nesbet in his report on Problem 234. 

219. COMPARISON OF SIMPLEX AND RELAXATION METHODS IN LINEAR PROGRAMMING 

A program has been written and checked for the solution of simultaneous equations by a method of » vein 
projections. If Ax < b, A an m x n matrix, and X is an approximate solution, * better solution, X t «i is 
obtained from 

v v r * e i mod m a T 

i+1 i . .T .2 t mod m 

(A e ) 

l mod m 

where e^ (05K frn-l) is the vector with (K+l) th component 1 and all other components 0, and 

r * max (r , 0), Ax - b « r, 
i ill 

In spite of the computational speed of this technique, the convergence is so slow as to render it imprac¬ 
tical. Various methods of accelerating the convergence have been proposed and are under consideration. 

A program for the solution of linear programming problems by a method similar to one proposed by 
Cooper, Charnts, and Henderson ( An Introduction to Linear Progra mming. J. Wiley and Sons, New York, 1954) 
has been written and is being checked. Other methods will be compared with this technique. 

This problem is being carried out by Mr. Dean N. Arden and Mr. Elliot Raiffa of the Digital Computer 
Laboratory Staff. 

221 B. COURSE 6.25, MIT 

Sixty-one students enrolled in course 6.25, MACHINE - AIDED ANALYSIS, given first semester, 1954, 
by the Electrical Engineering Department. The course was jointly given by Profs. W. K. Linvlll, R. C. Booton. 
and C. W. Adams. Computational techniques of analog and digital computers, the solution of engineering prob¬ 
lems by numerical methods, and numerical analysis were studied during this course. 

TAC, the summer session Three Address Computer, was used by the students for the solution of tne 
following problem: 

Consider the Van der Pol equations 
** 2 

y-» 6(y - 1) y ♦ y ■ 0. where 
y (0) « a 

y (0) » b 

For this equation, assumed values of £ . a, and b are chosen so that y < lb, y < lb, and y < lb. 

max max max 

Solve the equation using intervals, h » 0.1. Print out the values of y for integer values of t and carry oul the 
solution until the value of y has passed its second positive peak. Scale factor the solution so that all of TA< ' -• 
8 digits are used. Use the rectangular rule of integration. 

Each student programmed and prepared on punched tape his own program. Each solution took upproxl 
mately 3 minutes of WWI computer time and a total of 37 minutes was used by this problem during this report 
period. 
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222 C. HELICOPTER ROTOR STABILITY 

The Irons lent response in the flapping of a helicopter rotor blade has been calculated by Y. Shulman of 
the Mil Ac roe las tu and Structures Research Laboratory. These calculations involved the integration, using 
the fourth order Bunge Kutta method, of two ordinary, linear second-order differential equations with variable 
coefficients. One purpose of these calculations was to check the results obtained by an approximate analytical 
solution to this problem which was dealt with by Y. Shulman in his Master's Thesis*. The other purpose was 
to ' o m pa re results obtsined by assuming the blade to be flexible in bending to those for rigid blades and check 
with experimental results. 

Calculations were carried out for a variety of blade mass constants and over a wide range of forward 

speeds. 


The results have shown conclusively that within the limits of the approximations, blade flexibility affects 
the stability regions of the motion considerably. The results of the analysis which include blade flexibility 
tend towards better agreement with experimental results. The results of the analytical solution are uncon¬ 
servative. 

A full description of the problem, and the results of these calculations were included in a paper delivered 
by Y. Shulman at the 23rd annual meeting of the Inatltute of Aeronautical Sciences. January 27, 1955**. 

223 B. N. INVESTIGATION OF TURBULENT FLOW 

This project is sponsored by the Office of Naval Research under contract No. N5-ori-07874. It deals with 
the investigation of turbulent velocity fluctuations in open channel flow by means of a Pitot tube-pressure cell 
combination. In previous summary reports, this problem was described under problem 107. 

During the interim, the phase of this investigation concerning the measurement of some turbulence 
rharurteristu-B in the wake of a circular cylinder in supercritical open channel flow was completed. Auto¬ 
correlation curves were obtained from the Digital Computer Laboratory for points near the centerline of the 
wake .it three stations 40, 50. and 70 diameters downstream from the test body. From these curves, by the 
process explained in Summary Report No. 40 under problem 107, the scale of macro-turbulence or the mean 
eddy size was computed. The variation of this parameter with distance downstream gives an indication of the 
type of turbulent dc ay present in the wake. The values computed showed an increase with the square root of 
the distance downstream from the cylinder. This indicates that the rate at which the smaller eddies in the flow 
are decaying into heat must be greater than that at which the larger ones are breaking up. This same type of 
variation has been observed by investigators in the wake of a grid and is indicative of the isotropic nature of the 
turbulence present in the core of the wake. 

Additional runs are planned by F. Ralchlcn of the MIT Hydrodynamics Laboratory for the near future to 
ascertain the applicability of a piezoelectric ceramic transducer with impact tube attachment to the turbulence 
problem. 


224 N. COMPUTATIONS OF THE FIELDS OF VERTICAL VELOCITY AND HORIZONTAL DIVERGENCE 

ThU problem was prepared for compuUUon by Getrmundur Arnason and Is being programmed by William 
Wolf. The result* wtU b« analysed by the Preaaure Change Project under the supervision of Professor James M. 
Austin of the MIT Meteorology Department. 

The rising and sinking air motion causes the weather to change from clear Bkies to rain and vice versa. 
This important vertical motion cannot be measured directly so that the meteorologist must estimate the vertical 
velocity by observing Its effect on other measurable quantities, such as temperature. To date our information 
on the field of vertical motion Is very Inadequate. It is the objective of the proposed set of computations to 
determine this field for a series of weather situations. The results should contribute to basic tmowledge 


* Shulman. Y.. "Stability of a Flexible Helicopter Rotor Blade in Forward Flight" 
Engineering Department, Massachusetts Inatltute of Technology. June 19S4. 

** " aUblUt > of * n ” u >"- H.Ucop,er Rotor Hl.de In Forward Flight'. 


S. M. Thesis, Aeronautical 
preprint No. 521. IAS, 
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concerning the behavior of the atmosphere. Furthermore, they should greatly aid the eonUnuing research in its? 
field of numerical weather prediction. Various modified forms of Hie equation for the change in the vertical 
component of vortlclty are uttlMed to predict tomorrow's flow pattern. A lack of knowledge of the relative 
Importance of the vertical motion terms In the vortlclty equation is hampering the research in the field of 
numerical weather prediction. 

The main problem ia to determine the three-dimensional fields of vertical velocity (and horizontal 
divergence) by means of the observed pressure field. For this purpose the adiabatic equation is Integrated 
between the two pressures p • p t and p ■ p 2 (p ( >Pjl. 

log pj log p, 

(i) sT <z a' z i ,+ g J v -vt d log p * J wdiogp.o 

log P 2 log Pj 


Moreover the horizontal divergence ia determined from the continuity equation 


(2) 


div V • W + 


iW 

"h log p 


where Zj and Zj are the heights of the pressure surfaces p t and pj respectively. R is the gas constant for dry 
air. g the acceleration of gravity. V the geostrophlc wind. T the temperature, V T the horUontal nabla operator, 
S the static stability, and W the vertical velocity defined as 

,3, W - - ^20 

at 

d 2 

where — is the total time derivative and the local time derivative, 

at } t 

Observed Z's for p • 1000, 850. 700. 500, 300. 200, and 100 mb are available 12 hours apart in time. 

The first term in (1) is replaced by a 12-hour change in Z and interpreted as a local time derivative at the 
midpoint of the time interval. The geostrophlc wind V is related to the height Z by the equation 

(4) V • * K X VZ 

where f is the Coriolis parameter and K a unit vector directed upwards. 

Observed temperatures are not available, but an approximate temperature field can be derived from the 
heights (the Z-field). This is done in part 1 of the program. Since both V and 7T arc thus known, the second 
term in (1) can be evaluated. The third term in equation (1) contains the only unknown quantity W, since the 
static stability S Is derived from the Z-distribution (part 1 of the program). Nothing is known about the W- 
distributlon except that W * 0 at p ■ 1000 mb. An approximate vertical distribution of W is obtained by over¬ 
lapping second order polynomials in log EL of the form 

P 

(5) W • Wj ♦ log ( ^ }a(x,y) £log b(x,y) 


where a(x,y) and b(x,y) arc two parameters to be determined from observations, x and y are the horizontal 
coordinates. Wj is the vertical velocity at p » pj and W the vertical velocity at a variable p. We now put p * 
1000 and p 2 * 700 mb and carry out the integration of (1) utilizing the boundary condition W - 0 .it pj • 1000 mb. 
The resulting equation contains the two unknown n(x,y) and b(x.y). An additional equation is obtained by putting 
pj • 700 and p 2 » 500 mb and Integrating (1). We have now two equations with two unknowns and these are solved 
by conventional methods. Upon determination of a(x,y) and b(x,y) the vertical velocity at 700 mb is obtained from 
equation (5). Starting anew with p\ • 700 mb and W t - W 7QQ known, the procedure is repeated lp 2 . 500 and 300 
mb respectively) and a new set of parameters Is determined. Altogether four sets u, b are computed und the 
vertical velocity profile is composed of four overlapping parabolas with horizontal axis, see equation (5). 

The field of horizontal divergence is readily computed from equation (2). 
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Part II of the program deals with the computation of the quantities W and div V. The fields of vertical 
velocity and horizontal divergence being known besides the original Z - field, all the terms in the vorticity equation 

(6) ^ V -W — ♦ -I div V - PWx —— • K • 0 

3 t < i log p ' a log p 

b 2 2 

can be computed. Here ■ * V Z+f is the absolute vorticity, f the Coriolis parameter, and V the horizontal 

Laplace operator. Due to insufficient knowledge of the fields of vertical velocity and horizontal divergence, the 
three last terms have frequently been neglected in the applications of equation (6). 

The third and the last part of the program deals with the computations of the various terms in (6) and 
gives a simple statistical measure for the relative orders of magnitude. 

The programming of the problem is nearing completion. The program cycles through without alarms 
but not all of the answers are correct at present. It is expected that the conclusive checkout will be completed 
at the beginning of the next quarter. 

An interesting facet of tlie programming is the use of the symbol generator* as an output device. This 
special mode of scope output displays a coded combination of seven lines comprising a rectangular figure eight 
in approximately twice the time it takes to display a point. The usual type of number display being a coded 
combination of an 3 by 5 array of points, the use of the symbol generator effectively reduces the machine time 
by a factor of four. 

The coding employs the WW1 order code exclusively since the four significant digits required may be 
carried in a WW1 word. Where averages are performed a routine for combining the major and minor parts of a 
double length register number is used which is similar to that used in the Programmed Arithmetic routine. 

A general, yet concise subroutine for treating the Auxiliary Drum as a continuous storage medium was 
written. This was necessary since all of the available Auxiliary Drum storage (with the exception of about 200 
registers) is used during the course of the program. 

225 B.N. NEUTRON-DEUTERON SCATTERING. 

A report on this problem is given in section 2.2 of Part L 

228 N. EVALUATION OF DIFFERENCE DIFFUSION EQUATION 

A report on this problem is given in section 2.2 of Part I. 

230 C. DYNAMIC ANALYSIS OF BRIDGES 

of vi ° f 15 * MIT ^ lvU and SiniUry En * inecr ' n « Department 1. determining the dynamic re.ponae 

o variouH type, of .ample .pan bridge structure.. The purpose of the project is to determine the minimum value 
of a parameter q that is required to cause a predetermined maximum response. The parameter q defines the 
forcing function for any orientation of the impulsive load to the bridge. 

F ° r °' ,hU ^ nv estlgatl°n. a bridge ia represented dynamically by a Ringle degree-of-freedom 

system consisting of s concentrated mass and s spring, having various resistance .deflection characteristics. 

The solution is achieved by application of the second difference equation 

Ut) 1 *(t n ) ■ x(, ntl > - 2x(, n ) e x(t n 


to the differential equation 


m V it ) ■ P - R 
n n n 

where 

x(t ) ■ acceleration at time (t ) 

n n 

P • force applied at time (t ) 
n n 

R ■ resistance developed in spring at tilhc (t ) 

n n 

The overall program is advancing satisfactorily along the lines reported previously. 

231 B,N. REACTOR RUNAWAY PREVENTION 

A report on this problem is given in section 2.2 of Part L 

232 B. ENERGY LEVELS IN A SPHEROIDAL SQUARE WELL 

The problem being studied by K. Gottfried of the Laboratory for Nuclear Science is an investigation of 
Bohr and Matte Ison' a Model of Nuclei. A first step in such an analysis is the computation of wave functions 
and energy levels for nucleons moving in a deformed potential. The investigation of this first step 1 (which U 
described in Summary Report No. 40) is virtually completed, and a final report will be submitted in the next 
Quarterly. 

233 C. UTILITY STOCK PRICES 

This project, which involves statistical regression analysis of utility stock prices, is being conducted 
by D. Durand of the School of Industrial Management with two distinct goals in mind. In the first place, it is an 
experiment in applying regression analysis to financial problems in order to (1) develop better methods of 
ascertaining the cost of equity capital and (2) to establish principles applicable to corporate financial manage¬ 
ment and to security analysis. In the second place, it is an experiment in reducing the arduous computing labor 
heretofore encountered in regression analysis to practical proportions by the specific means of a modern, tape- 
fed electronic computer. 

The problem entails a series of computations that can be classified into four groups, as follows: 

1. Preliminary processing of data. 

2. Summntion of squares and products - an<1 

3. Solution of simultaneous linear equations with concurrent operations required 
in statistical regression analysis. 

4. Final processing. 

Of these, the nature of 1. and 4. is not yet decided. Possibly all or port of the operations required can be per 
formed by hand or with desk calculators. One preliminary operation that is not well suited for hand work 
is the conversion of some five to ten thousand decimal numbers Into logarithms. Group 2 ia laborious for desk 
calculation but conceptually simple; a subroutine for doing this work on WW1 may easily be set up. Group ' 
presents the real challenge. 

Although the literature on solving simultaneous equations is voluminous, there appears to be room for 
considerable development work. In statistical regression, it is necessary to solve the so-called normal equa¬ 
tions. which could be done by any of the already developed methods, but the problem does not end here. In 
addition, residual variances and standard errors have to be calculated. Although these- additional operations 


•Por a detailed description ace M-1623-2 "Programming for !r. Cut Unita". 


'See also K. Gottfried and V. F. Welaakopf, (Quarterly Progreaa Report, L.N.S.. February 1#55 








APPROVED FOR 


WHIRLWIND CODING AND APPLICATIONS 


. in be treated separately after the solution of the normal equations, they can also be incorporated directly into 
the solution of the normal equations. The scientific contribution, if any, will come from the effective integration 
of a number of computation steps into a single systematic routine. 

Progress to date has consisted in writing an experimental program, which is now in the process of being 
tested. At the same time we arc thinking about possible improvements in this program that will make It more 
grnrral and extend its usefulness to a wider class of problems. 

234 N. ATOMIC INTEGRALS 

Electronic wave functions for atoms and molecules (for states of arbitrary symmetry) can be obtained 
to arbitrary accuracy by a method which consists of three relatively independent stages of calculation. 

Stage A - choice of some finite set of single-electron functions, . and evaluation of all possible one- 

and two electron integrals over this set. The integrals are matrix elements of operators occurring in the 
many-electron Hamiltonian. 

Stage B calculation of expansion coefficients of a set of orthonormal single-electron functions . 

Stage C calculation of configuration interaction effects and resolution of degeneracies in the variational 
matrix for the many - electron Hamiltonian. 

If Stage B is carried out. with certain modifications which are too complicated to explain here. Stage C 
is greatly simplified. Techniques of group representation theory and perturbation methods can be used. The 
actual amount of calculation required is very small compared with Stages A and B. 

The eMsenUai difference between atomic and molecular calculation is that in the atomic case there exists 
a set of basic functions . the analytic Slater orbitals, which lead to rapid convergence in expansion of the 
mo If-conflate nt • orbitals, and for which all integrals can be evaluated in closed form. No class of functions is 
known which has both these properties in the molecular case. 

Programs are available at present which carry out Stages A and B for atomic wave function. These are 
not yet in their mo.st efficient form and some further programming is needed to join them into a single unit. 

Dr R. K. Neabet of th»* Solid State and Molecular Theory Group has programmed the calculation of integrals. 
Stage A for atoms and the program for transforming integrals under Stage B. These programs are described in 
detail elsewhere. 1 


orbitals 


The program for atomic Integrals re art. In a list of parameters ( A, Z^. specifying a act of Slater 


A + i 


-Z r Y‘" (8, g). 


oner P rr? °“' ‘"? Cp ' nd * W 0 ™'* nd ‘^electron Integrals (overlap, kinetic energy, and Coulomb potential 
onerg>). for normalized *\ a. The normalization constants or- also printed out. Here A and / are anv non 
negative Integers and is any positive number. • y non 

problems!* l“ie7,roT bleTaf ^nTk.T? B “* Uy lhc ■» both atomic and molecular 

1 and J or of k and I or of the indi ps'Juj, X I '«*£“T teK “ C ** 

** i 


*•« ••P"' 1 

the transformation required for two-electron integrals Is 


‘" * *>»■ ■..» ... 7*.„„ Group MIT. April 
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H*']- * f f f x -. x ,j x »k x „ L»i k *]. 

All functions and coefficients are assumed to be real. The corresponding transformation on two Indices for one- 
«• U*i iron integrals is considerably more manageable. 

The Whirlwind program for this transformation constructs all independent elements from 

[i j | k IJ without repeating equivalent calculations. The symmetries of this four-index form under interchange 
of the indices as well as any symmetry properties of the basic orbitals (due to spatial symmetry of the Hamil¬ 
tonian) are taken into account to reduce the number of arithmetic operations in this transformation as much as 
possible. 

23* B,N. EIGENVALUES FOR A SPHEROIDAL SQUARE WELL 

A report on this problem is given in section 2.2 of Part I. 

236 C. TRANSIENT RESPONSE OF AIRCRAFT STRUCTURES TO AERODYNAMIC HEATING 

The study of transient response of aircraft structures to aerodynamic heating was initiated on January 3, 1055 
by L. A. Schmit of the MIT Aero-elastic and Structures Research Laboratory. H. Parcchanian has been responsi¬ 
ble for a major portion of the programming during the latter half of this report period. 

The overall problem is that of investigating the influence of aerodynamic heating on the structural design 
of high speed aircraft. One important step in the solution of the overall problem is to determine the transient 
temperature distributions in built-up aircraft structure. This first phase requires the solution in generalised form 
of two idealized heat flow problems. The first of these two idealized heat flow problems (Problem I) has been 
formulated and solved. 

Problem I is that of determining the transient temperature response of a thin plate exposed to a timewise 
step function change in adiabatic wall temperature with a turbulent boundary layer type chordwlsc variation of 
the heat transfer coefficient. 



FIGURE l 

SCHEMATIC REPRESENTATION OF PROBLEM I 


The problem Is formulated in terms of the following symbols. 


T 

■ adiabatic wall temperature 

aw 


Ms) 

• heat transfer coefficient at x 

JL 

• plate length 

X 

• plate thickness 
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h 

0 

X 


heut transfer coefficient at x » 0 

local iuii along plate referred to leading edge of plate considered 

(5) 

V 

ait f + r 

Mi) 2 L 

1 f.lilai - 2 A ‘£1 f. 

' (fix. 2 J ‘•J- 1 

(Ax) 2 w ')' 1 

t 

T 


time 

Temperature 

(6) 

*n.j ‘ 

Mi) 2 W« 

i-r <x)At - « 

" Mx) 2 n ‘J- 


/ 


weight density of plate material 






k 


thermal conductivity of plate material 

or in matrix form 





C 

P 


specific heat of plate material 

(7) 

N 

■ w h-.l 




■ thermal diffusivity of plate material 


k 

PC 
' P 

If it is assumed that h(x) is nut n function of t or T, and the thermal properties of the plate are constant, the 
heat balance equation for the plate shown in Figure 1 is as follows: 


II) 


h(x) (T 


'** . dx 2 J ~P ' « 

The intttial and boundary • onditions are assumed to be as follows: 
(2) T « T k when t - 0 


(3) 


dT 

dx 


0 ut x = 0 and x • i 

The problem is put in nondlmcnsional form by introducing the following additional notation 

h 


T T 
aw _ 

T - T. 
aw i 


t « 


/■ s 


and 


r<s> • ~ 

n 


fin * 


With the notation above, the problem is stated in nondlmcnsional form as follows: 


(la) 

d# 



— + 

f(x)0 * fi 


dt 


(2a) 

0 * 1 

when t * 

13a) 

dj 

dx 

0 when x 


d 2 » 

di 2 


i-.i- - ? - . 

formulation can be written as follows ’ * he ° uncoupled forward finite difference 


( 4 ) 


'll * At A—*] » e Alt. , 

• J L 1 <a *> 7J i.i-t "j.j-t 


with the initial condition 


(<.)■ H 


The final production program is based on a fifteen element physical grid (n*15). Given numerical values 
for the parameters ^ and f(l) as well as a value offlf the program computes the coefficients that make up the 
matrix £Cj in eq. 7 using CS 11. Since it can be seen from the formulation that all non-zero terms in the matrix 
[c] are between zero and one the results of the CS 11 calculations are stored as single length WW1 numbers. The 
evaluation of |0j] is then carried out according to eq. 7 in WW1. Since the bulk nf the calculation consists of per 
forming the matrix multiplication indicated in eq. 7 the time saving that results from performing this part of the 
program in WWI rather than CS (1 is substantial (at least four times faster). 

The output display is controlled by the nondimensional temperature response (1 0| j) of the leading edge 
element. As (1 -01,j) reaches 0.1, 0.2, . . . etc. up to 0.9 at intervals of 0.1. the fifteen element nondimensional 
temperature distribution is stored as well as the time (t jAt) associated with it. When (l-0j j) exceeds 0.9 the 
stored results for the case <f( 1 >. just completed are displayed. On the first scope frame. 1 0i # j versus x for 
various t are plotted. The first frame is essentially a family of nine nondimensional temperature distribution 
• urves. On the second scope frame 1-0 jj versus t for various x are plotted. The second frame is essentially 
a family of fifteen nondimensional temperature time history curves. The nine pertinent times are typed out via 
magnetic tape. The final values of 0,j are also typed out via magnetic tape and would serve as an initial condition 
in the event that it ever became necessary to carry the calculation nearer to the steady slate condition. 

The selection of At in order to prevent divergent oscillation of the difference solution is a routine matter. 
(Kef. F. B. Hildebrand. Methods of Applied Mathematics, Prentice Hall, Inc., New York, 1952, pp. 323-345.) The 
maximum value of 211 that may be used for a given p in the fifteen element program is given by eq. 8. 


( 8 ) 


A t 


1 


max 450 p t 1 


It was found by trial that satisfactory damping of the convergent oscillations, characteristic of the finite difference 
solution at early times, could be insured by selecting A t such that approximately four hundred time cycles are 
required for the nondimensional temperature of the leading edge element (l-0i j) to reach 0,0 of its steady state 
value which is 1,0. 

The convergence of the fifteen clement solution was examined by writing a forty-five element program and 
running cases at the limits of the ^ and f(l) parameter ranges. It was found that the maximum difference between 
the fifteen and forty-five element response occurred when f(l) was a minimum and p was a maximum. (f(l) ■ 
0.1, £ • 5 x 10' 2 ). As a result of the calculations made employing the forty-five element program it was decided 
that the fifteen element program yielded results satisfactory for the purposes of the present work. 

During the course of the work Just discussed it became apparent that for certain values of ft and f(1) 
satisfactory results could be obtained by neglecting chordwlse conduction. The nondimensional temperature 
response is obtained tn simple closed form if chordwlse conduction is neglected. 


(9) 


d0 


4 f(x)0 * 0, 0 ■ 1 when t ■ 0 
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ftecuuso of thu fact the program has been augmented. For each case (f(l),P) there are 135 values oil a 
storage Just prior to entering the output display portion of the program. The augmented program computes the 
136 corresponding values of when conduction Is neglected, according to eq. 10, then computes and prints 
out via magnetic tape 135 values of e, whore e is determined as follows: 

t - ft .1 conduction neglected 

(11) e-JiJ - ‘ 4 - - - 

' Vi 

As a result of running 33 cases employing the augmented program the following tentative conclusion has been 
reached. If ^ < 5 a 10 '2 jpg 0.1 * f(i) < l.o, less than 5% difference In the nondlmenslonal temperature reBpons. 
results from neglecting chordwlse conduction provided 

(13) fill > 0.45716 log ft + 1.43347 

Work on Problem I Is essentially complete. It is intended to program Problem II. the web-plate tempera¬ 
ture reBpons.- problem, during the nest report period. Formulation and discussion of the results obtained for 
Problom II will be presented In the nest report. 

337 C. AUTOCORRELATION FUNCTION OF SUBMITTED DATA 

This problem was concerned with finding the autocorrelation function of data derived from work originating 
at the Instrumentation Laboratory at M.l.T. The problem used the subroutine developed by Douglas Ross for 
computing autocorrelations In problem 107. The autocorrelation function was derived and the problem was 
completed. 

338 B,N. SELF CONSISTENT CALCULATION OF NUCLEAR MASS DENSITY 

A report on this problem is given in section 2.3 of Part 1. 

239 C. GUIDANCE AND CONTROL 

Th ‘“ p ™ blem encompasses a variety of calculations relating to the guidance and control of aircraft. De 
tails are for the most part classified. However, as incidental by-products of this investigation, standard sub- 
routine library tapes have been prepared for: " UD 

1 ) solution of an arbitrary number of ordinary first order differential equations using the method of 

2 ) the minimisation of an arbitrary function of n variables by the method of steepest descents. 

Further discussion of these subroutines may be found under problem 141 . 

tion LabLrato^ W,1i * atl0n “ h'*"* ° U ‘ “ nder ,he ‘^"islon of Dr. J. H. Laning. Jr., of the In.trum.nla- 

241 B,N. TRANSIENTS IN DISTILLATION COLUMNS 

members «»rk has been used by 

too complicated to solve by ordinary methods Up to the'nres.V u problem ® <H««Hlatlon which are 
solution of some complicated 1 J W ° r b<?e " d ° ne by Jack °’ Donne11 on ,hl ' 

been don. b, Smith. Polk, and Jordan on unsteady state ^ ^ *** 
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All of the work on unsteady-state problems has been on an ideal system which satisfies the following 
conditions: 

1) A two-component system. 

2) Constant molal overflow above and below the feed plate. 

3) Equilibrium condition between the liquid and vapor phases represented by y • — ^ where x and 

y are the mole fractions of the more volatile component in the liquid and vapor phases respectively. 
Furthermore at is assumed to be constant throughout the column. 

4) Average composition of the liquid leaving each plate equal to the composition of the liquid on the plate. 

5) Average composition of the vapor leaving each plate given by the equilibrium condition with the liquid 
on the plate. 

6) Holdup of all the plates constant with time and equal for all plates. 

7) Holdup of condenser and reboiler negligible. 

8) Feed introduced into column as a saturated liquid. 

The operation of such a system can be described by u set of non-linear, first-order differential equations. 
Such a set of equations can be solved by any of a number of different types of finite-difference approximations 
to almost any desired accuracy. The Runge-Kutta approximations have advantages which make them the easiest 
to use with a large scale digital computer. It has been found that by uBing a second-order Runge-Kutta approxi¬ 
mation with an interval aize equal to 0.1 units of liquid throughput (time x liquid flowrate /individual plate holdup) 
insures convergence and at least 0.05% accuracy for nearly all problems of this type. 

Programs have been written for the CS II computer for a distillation column containing up to fifty plates 
in which any one or more of the following variables have undergone a sudden change: 

1) Feed composition 

2) Reflux ratio 

3) Vapor rate 

It is proposed to study how a column will react to such sudden changes and to attempt correlations for the 
length of time that it takes a column to approach its new equilibrium condition. It is also proposed to write pro¬ 
grams for and to study a column in which some variable is controlled to offset the effects of a change in some 
other variable. For example the reflux ratio might be controlled while the feed composition undergoes u change 
so the tops composition remains approximately constant. 

Finally it is proposed to write programs for and to study more complicated unsteady-state distillation 
problems in which one or more of the above eight conditions are relaxed. 

These studies will be carried out by S. H. Davis. Jr. of the Chemical Engineering Department. 

242 N. NUMBER OF STRUCTURES OF RELATIONS ON FINITE SET 

A report on this problem is given in section 2.2 of Part 1. 

243 D. CRYSTAL FILTERS 

The data obtained from this problem is designed to provide a simple approximate solution to a phase of the 
design of filters using quartz crystals as elements. The part done on WW1 consisted of the systematic computation 
of products, using cycle counters, with the results being obtained on film. There are no plans for nny further work 
to be done on WW1. The problem was done for David Kosowsky and will be published in his Ph.D. thesis and in the 
Research Laboratory for Electronics Technical Report. Hannah Paul was the programme**. 

244 C. DATA REDUCTION FOR X-l FIRE CONTROL 

This problem is concerned with computing from fire control signals how far fictitious projectiles miss an 
observed target. Because of the large amount of data taken for each target run. use is made of the great storage 
capacity of the auxiliary drum. The main tape for this problem calls data off the drum and Into high speed storage, 
a block at a time, for each cycle of computation. 
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Calculation* when- accuracy ia a factor, in general, ore carried out with programmed arithmetic where 
aa Home .cac hing procedure, are coded in the basic code of WWI. Extensive ballistic data also is stored on the 
auxiliary drum. Register numbers for needed ballistic data are computed in the basic code of WWI, 

A nu ' nber of 8 Sorter tapes prepared for this problem are concerned with analyzing system errors and 
finding suitable corrections to apply to dial readings for the raw data used with the main tape. 

This work is being carried out by Dr. J. M. Stark of the M.LT. Instrumentation Laboratory. 

245 N. THEORY OF NEUTRON REACTIONS 

Recent treatments of the theory of nuclear reactions have shown that as far as neutron reactions are ™ 

^tentui taVta f’l Tr'V' ,he P / 0bl,;m COUld b * explam ' d b > assuming that the neutron moves m a complex 
CSTrWtatodt first treatment of this problem it was assumed that the potential was in the form of a recUngular 
ell. ThU led to relatively good agreement with the total cross-section* and to only rough agreement with th 

r«:«: zT rlns * nd r uUr du,r,bu,ion ° r 

weUm'^ '" tOUrU8<! , d “ l ° °°- “ WaS fe “ lh “> 8 Pufontial w.U which was no. as idealized a.t^snulre' 
weu. in particular a weU which had a tail ,o «. would bring the experiment and the theory into cl^r axemen,. 

de.crip£n llm ‘ n * ry '*° r, ‘ *" md ‘ CaUd ‘ *' U «* form should be close to the correct 


V ■ -V o (W,$> [i . unh^-l-ij. 


The parameters describing the well are 

V D * the depth of the well 
R * the approximate radius of the nucleus 
d -the rate at which the square well falls off to zero 
^ » a measure of absorption inside the nucleus 

we,, of this Zt'S .Tm U uU,:l en r::r.Te d ,^rrmel r Z\°s “ » h “*> “ - «— that a 

into . nondimenslonal form by uitroducing the m ' V Pr0t ° n8 - The problem *« transformed 

« 2 2m „ _l . 
o £2 V o R d ^ R » and x ■ kR 

-here m ■ m... of neutron and k 2 . 2m/* 2 E where E is th. neutron energy. 

the best value, of / and { so a. to LstTpe^m/rt wit”£u. °* ViaUe ’ °' P ‘‘ r * n ’* ter " «° P‘ ck out 

The radial equation to be solved for the J < b partiai wave U 




I Y 


u^ (0) . 0 


where the potential 
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A 

V(y).l/a—j- (1+iCl (1-tonh 

From U|, the phase shift jjj can be found by the equation 


J y n(y> V<u) u. dy - 1.3 . . (2 Jsl) 

«* { ■ --- r- -Z2_zao 

J y Jjj (y) V(y) Ujf dy 


From cot can be found 


% 


/ rot 2 / -1 cot J. 

.e*4._i_ + 21-=--- 


cot 2 sj 41 cot 2 «£ 41 


will give all the J th partial cross-sections (total, reaction, and elastic scattering). 


^*^ 2 V V 


rr r 


•n R* 


* aJt 

_2 


,J r a -n, 1*1 


21-tl 

2 


I 1 M* 


A program haa been written which finds the value of the wave func tion Mg ■ Vj -t twj by mcona of a power 
aeries for w and v. Another program is being worked on which will generate the wave function from these start¬ 
ing values by using the Milne Predictor-Corrector formula. At the same time thr function is generated the inte¬ 
grals in the expression for cot will be computed. As soon as this program has been tested it will be combined 
with the first one in such a way that given any x. X Q , S and / the cross sections will be found for A from 0 to 6. 

This problem is being programmed for Professor H. Feshbach of the Physics Department by E. Campbell 
and E. Mack of the Joint Computing Group. 

247 C. SURFACE PRESSURE PREDICTION 

In the past a considerable amount of work has been done by Prof. Wadsworth’s D.I.C. Statistical Laboratory 
at M.I.T. in the prediction of surface pressure by means of linear functions of present and past values of pres¬ 
sure, taken over a network of points on the map. As a collateral study by Prof. Wadsworth's group, tt was found 
that a large sector of the pressure map could be represented adequately us a geometrical surface, expressable 
analytically as a general polynomial In the coordinates of the points, regarded as lying on a plane surface. By 
this device, it was possible to reduce a network of 91 observation points to a set of 15 orthogonal function*. 
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namely, the • ocfflcu*nt« of the polynomial expression. Using these orthogonal functions in place of direct pres¬ 
sure observations, Prof. Malone carried Wadsworth' * method of pressure predu tlon a step further, with very 
em ouritging results (cf. problem 155). 

An entirely independent approach to this problem has been made by the so-called numerical forecasting 
group. The technique employed by this group involves the numeric al solution of differential equations with stated 
boundary conditions. The type of solution depends upon the physical assumptions made in setting up the differen 
Dal equations, but in general the form of solution in largely linear in the pressure values. However it does con 
lain . rUin quadratic terms. The results of Prof. Malone proved to be at least as good as those obtained by the 
numerical forecasting group, but it seemed desirable to exploit the idea of quadratic functions further. 

The object of the present study by the D.l.C. Statistical Laboratory is to evaluate the introduc tion of the 

• tuadratn component into th« general linear scheme of Wadsworth and Malone. The pressure surface is repre¬ 
sented by .* et of orthogonal functions, and then the squares and products of the most important of these functions 
are incorporated us representative of the general quadratic contributions. 

249 T. FLIGHT INTERCEPTOR CONTROL 

This problem is concerned with a systems-analysis investigation of the dynamics and control of an inter- 

• eptor. The system is described by thirteen differential equations, which were solved on Whirlwind l by the GUI 
method. The Whirlwind 1 solution has been completed, and its digital results served as a standard check for later 
analogue * omput.itions done on the M.I.T. Flight Simulator of the Dynamic Analysis and Control Laboratory. 


By convention, the positive direction on a straight line is specified by the order of the two points uBed 
in tagging it. and is the direc tion from the first point to the second point. The positive direction on a circle is 
always counter clockwise. Starting from the point of minimum x and proceeding along a curve in its positive 
direction, the first intersection with a second curve is the near IN) intersection; the other is the far (F) inter¬ 
section. Two circles, or a line and a circle are said to be tangent (T) when they have one and only one inter¬ 
section and their positive directions coincide at the point of intersection, they are antitangent (A) if their positive 
directions do no coincide at their only point of intersection. 


Typed Symbols 

P 16 • -2.734, 6.2545 
p35 - si2 |sl7 
pl2 • N si | cl2 ] 

pi2 = F al| el5 J 

plS » N c3| cS 1 
pl5 * F c31 cl6j 
pi9 = ell | 82° 


coordinates 

intersec tion of two lines 
intersection of a line with a circle 
intersection of two circles 
on circle, angle with positive x axis 


250 1 . TRANSLATION PROGRAM FOR THE NUMERICALLY CONTROLLED MILLING MACHINE 

The M.LT. Numerically-Controlled Milling Machine (NCMM) is a standard milling machine modified so 
that it i an be automatically c ontrolled by numerical instructions punched in a paper tape. A comDlicated sequence 
of motions of the cutting tool may be entirely specified in speed and direction on the punched tape, and no manual 
intervention by a human operator is required during the cutting. 

The data punched in the control tape must appear in numerical form. This requires, in general, that the 
pret ise x, y. and z loordmates of all significant points on the work be computed. The tape instructions specify 
the motion of the .enter of the cutting tool, while cutting occurs at the periphery of the cutter and not at its 
• enter. This necessitates a further transformation of the points on the work to account for the offset of the cutter 
' enter. These calculations, even for a relatively simple cutting job, are tedious and time-consuming. The task 
of preparing a control tape is even further complicated by the fact that the tapes must be punched in an octal 
code in which the final instructions take a form quite different from the numbers used during their calculation. 

In order to reduce the time required for the preparation of NCMM tapes and to minimize tape errors, a 
procedure for utilizing a digital computer has been found feasible. A translation program is being written for 
Whirlwind 1 whir h accepts a desc ription of the* work in symbolic form and which produces automatically the 
required NCMM control tape. The calculations and transformations which were described in the preceding 
paragraph are executed automatically by the Whirlwind I computer. 

The first form of the translation program has, for simplicity, been limited to curves consisting solely 
of straight lines and circles in two dimensions. Straight-line motion in the third dimension has been Included, 
in a somewhat restricted form. 

A brief description of the vocabulary accepted by this program follows. All symbols are typed on the 
M.I.T. Flexowriter exactly as they appear in the description. 


s3 ■ pi I p2 
s2 ■ pi | Tc2*| 
s2 * pi j Ac 3 J 
si = Acl | Ac3 
s2 » Tel j Ac3 ( 
s3 ■ Acl | Tc3 
s4 * Tel | Tc3J 

s3 = pi 2 I 73° 


c7 • pi. 2.7405 
c2 » pi | Tell 
c2 * pi | Ac3j 
c5 * p5 | s2 

Cutting instructions 

After the significant points, lines, and circles on the 
be specified in terms of them. The cutting instructions take 


two points 

through point, tangent to circle 


tangent to two circles 


point, angle with positive x axis 


center, radius 


center, tangent to another circle 
center, tangent to line 


work have been tagged, the actual cuts required will 
the following form: 


Points, straight lines, and circles will be tagged using a tag assignment of one of the forms given in the 
table below. A tag* may not appear after the * symbol unless it has already been assigned. 


•A tag consists or the letter p(polnt). efeir. le). or straight line), followed by any integer 1 such that 0 < 1 <255. 
1 he integers serve only to distinguish one tag Irom another, and have no numerical significance. 


.These Ug ' assignment, may be typed a. .2 * Tc2 pi and »2 • Ac 3 |pl in order to reverse the posits direction 
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Symbols Typed 

Meaning 

7.5. p3, -8.214 

move at 7.5 in/mln feedrate in a straight 
line to p3, at the same time lowering the 
cutter 6.214 in. 

15. tel, p2. 0.3 

move at 15 in/min along circle cl in its 
positive direction, to p2. at the same Ume 
raising the cutter 0.3 in. 

7.5, -c3, p5 

move at 7.5 in /min along circle c3 in its 
negative direction to p5 without vertical 
motion of the cutter 

STOP 

stop the NCMM (punch feedout) 

END 

signifies end of Flexowriter instruction 
tape 


The Initial cutting instruction assumes the tool is at x • 0. y 
to be where the preceding Instruction left it. 


0. Subsequent cutting instructions assume the tool 


Special Words 


Special words, used to 
on the input tape as follows: 


provide information needed in calculating the 


NCMM instructions, must appear 


RIGHT 

LEFT 

TOOL RADIUS - 0.5 
TOLERANCE . 0.0005 


Tool to right of cut 
Tool to left of cut 


‘STS?’ “ “ “^" d «•«* «. routine wlU be 

the Digital Con.pu.er Ub^rlry «.« ^ ^ ^ ^ •*"* U b * 1 "* —led out by A. Siegel of 

251 B. DYNAMICS AND CONTROL OF PACKED DISTILLATION COLUMNS 

of a fairly large number of phy^ral’para^Aera" m ” ° rd ' r to dc " rml,,e which, out 

dynamics of packed column* affecting ,hc composition at the .o„ r ^ USItd *° d " cribe and P redtct ‘he 

grammed to aimulale IhTjrtml Mrto!i <> ! ^"?I.* < |U * tt ° n * WUch h * ve be *' n dev ' 1 °P' rl - Whirlwind I la being pro- 
Ume.* e# * e d ^P'PP^^kruTctectaB^aga^iflt^measured^data.'T^Ua^^LM tah^n^^Qyj two^hour^of^computer 


The difference equation, are concerned with the following phyMc.l effect.: 

1) mass transfer; 

21 heat transfer, 

3) fluid and vapor flow. 


The complete program has not yet been completed, and therefore no conclusive results can be stated al this 
time. The study is being carried out by H. M. Teager of the Electrical Engineering Department. 

252 N. ANALYSIS OF TWO STORY STEEL FRAME BUILDING 

A report on this problem is given in section 2.2 of Part L 
256 C. WWI-ERA 1103 TRANSLATION PROGRAM 

A two pass input translation program has been written by members of the Digital Computer Laboratory at 
M.LT. to simplify the machine coding problem for programmers working on Problem 126. It will allow them to 
code programs for the Umvar Scientific Model 1103 in a mnemonic (two letter) operation code and to use either 
or both symbolic and integer addresses. The translation program accepts Flexo-coded punched paper tape as 
Input and produces a seventh-level bi-octal punched paper tape acceptable as input to the 1103 computer. The 
program operates on the Whirlwind 1 computer at M.I.T. and was written under the sponsorship of the Digital 
Computer Laboratory and Project DIC 7138 of the Servomechanisms Laboratory at M.I.T. It will be used in 
the coding of the large data reduction program which is being developed .it the Servomechanisms Laboratory 
under Problem 126 and which is expected to be run on the Egiln Field 1103. 

The 1103 computer is a large two address computer containing 1024 or 4096 registers of random access 
memory and 16384 registers of directly addressable drum memory. Each register consists of 36 binary digits 
of which, in most instructions, 6 bits are used to designate the operation and two sets of 15 bits designate the 
two addresses. In some cases the two addresses can be divided so that as many as four "addresses" can occur 
in the instruction. The registers in the arithmetic element of the computer are directly addressable. Both 7- 
hole punched paper tape and punched cards can be used for input and output. 

The following is a brief summary of the vocabulary and syntax of the translation program: 

A. Characters 


The words in the vocabulary of the translation program will be composed of suitably punctuated and 
terminated syllables of letters and digits punched in standard Flexowriter code on paper tape. 

B. Words 

Three classes of words are defined. The first class consists of all those words which occupy storage 
registers in the translated program. These words are instructions, with zero to four addresses, or (integer) 
numbers. We call this class of words polysyllabic storage words. 

The second class, called polysyllabic control words, consists of current address assignments, symbolic 
address assignments, and starting address assignments. These words influence the form, location, and operation 
of the translated program but do not themselves occupy registers of storage. 

The third class of words, called special words, are used for miscellaneous control and identifying pur¬ 
poses. They are the title, the number base indicator, and the comment word. 

The words of the first two classes are combinations of these syllables: operations, symbolic address tage, 
integers and the literal addresses "q," "a," and "b." These words are distinguished solely by the Internal and 
terminal punctuation. These characters are 'V', *V\ " | ", the tab. and the carriage return. 

C. Syllables 

1. Operations. These are the mnemonic lower case two-letter pairs corresponding to the standard 
listing of the 1103 operation code. 


" * - - - —r XZZSXEJSZZ 


numerical model involves 
of both time and distance 
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2. Symbolic address tags . Theac arc three character digit and letter combinations. 

3. Literal addre sses . The letters "q". "a", and "b" refer to the quotient, right accumulator and left 
accumulator, respectively. 
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4. I nt e g er., Integer, may have any value up to 2 3s -l, ao long as the completely evaluated «ord or o,r. 
word ' tnte ‘' r “ “ 3yUab, °- '“ 8 ‘ n, ° ,h * numb ' r of biu tneanmgtully available to the trammed 

D. Polysyllabic Storage Worda 

1. I nstruction. . An Instruction word always ha. the operaUon code as the Initial avllable Zero 
wordTl"*’ m “ y Ih ° ,r,n ‘ UUon P ro *ram will decide where to position the addre.se. blt-wue In the Ur 

termited - addr ° B “' 8 “ ^ •<*«.„. are ..plated by comm^ and thV U.t t. 

address., and lm °' ?*" ** " ,Urn - E,ch “ ddre8a “ wrl,te " “» ‘he sum of symbolic address tag. literal 
by a plus or rntHrSn m, * nner °' 8umm ‘“ lon u u ' d ‘ c '“« d by prefixing each syllable of the address 

- “ —— * - 

E - Polysyllabic Control Words 

tlon word except that here I’‘p* ^ r ° rm d ° f uddre “ s of an mstruc 
registers of storage untU another current a^Ta. La^t" WU ‘ ,hereaf *' r *° ln, ° 

*. s&smz .—-- - 

If one were the next word on the Flexo tape tW. ,s true whether lr „ , y ‘ y " OUW ° CCU P> 

several times In a program, the last valJTiwing usedlTte I„sTamd ° rCUr ’ th ' r '- A »*« be reassigned 
caled as such by the translation program. * translated program. AU reassigned tags are indl- 

Flexowntor u p T; FsJtS^rm rfTlmglfu "° rd "l*"* OCCUr at thc end of each 

-—- - smarts , szz,zr ,m 

f • Special Worda 

L A ,i,le mu «* occur at the beginning oTeluI'” UP *“ ^ l0gBlng th * tr " nslatl °n of the tapes on Whirlwind 

“o^ ta .m rc n o^'^Cin Wher ? “ “ “ ny * n,e «" f ™» «» to ,e„. 
indicator occurs. A "base 10" la assumed initially C ° nVer,ed *° “**"7 'com this base. until another base ' 

JtoJtowSlT^i °" “ pro ‘' ra ^p”I»hI'^ sod i. provided so that com- 

“ y t ' rmlW ‘ , ' d * ■ “ ° r —- return. Any’ Flexo 

G. Symbolic Addreaaea 

If a programmer decides that he — J4# 

f;f=“=~sfF»=s~=- 

r .^ & —^ :sirr.r.r. : —~r 



WHIRLWIND CODING AND APPLICATIONS 


programmer haa a complete record of the storage locations used by his program, as well aa of the effect of 
modifications on the program. 

H. Error Detection 

Illegal combinations of characters are detected, and so recorded, by the translation program wherever 
possible. The locations of the errors are given in terms of the most recent symbolic address or current address 
assignment. 

I. Present Status 


The basic elements of the translation program are now working. Translated program tapes are complete 
except for visual titles, but the programs for printing the symbolic address tables and the detected programmer's 
mistakes have not yet been tested. The WW1 utility input program must also be modified In order to allow the 
translation program to operate in the aamc automatic manner as CS and the other WW1 systems. All of this work 
will be completed In the next quarter. 

258 C. DYNAMIC ANALYSIS OF AN AIRCRAFT INTERCEPTOR 

This problem involves systems analysis of an aircraft interceptor control system and is being conducted 
to determine the adequacy of certain force and moment equations to describe the motion of an aircraft. The 
system is described by eight differential equations to be solved by thc Gill method. Four solutions using dif*' 
ferent types of inputs to the control will be made by Whirlwind 1. These solutions also will serve as a standard 
for analogue computations to be performed on the M.l.T. Flight Simulator of the Dynamic Analysis and Control 
Laboratory. 

The general program for the first Whirlwind solution has been written and tested and is now ready to be 
run. The second and third solutions will use the same basic program with minor changes necessitated by the 
changes of input. The fourth solution is being programmed separately and is now ready for testing. 

280 N. ENERGY LEVELS OF DIATOMIC HYDRIDES 

A study of the electronic energy of the OH molecule including configuration interaction is under way. The 
energy will be calculated as a function of intcmuclear distance. The ground state of the molecule is 2 7T. 

The basic set of one-electron orbitals are the hydrogen Is ground state function and the is. 2s, and sp 
Hartree Fock 1 wave functions of oxygen. The notation followed Is to call oxygen Is ■ s. oxygen 2s « d\ oxygen 
2 p * p, and hydrogen Is * h. 

Use is being made of the Whirlwind I digital computer to calculate one- .and two-electron Integrals of the 


form: 


“|*ih> ■/rV Uf i f j (Udr » 


,«V U ^» dl , dr 2 


“j |g 12 |x»■/**,«> f,' 

Uj|i 12 | xi * g l2 f t <i) f k <2> dfj dr 2 


f * b are the one-electron wave functions and fj and g l2 the one-and two-electron operators respectively, 
grain are easily calculated by using the program written by F. J. Corbato? The most difficult of the 


where the^' 

These integrals L _ 

integrals takes about 5 seconds of machine time while the others take considerably leas. 


l *D. R. Hartree, W. Hartree and B. Swl.les, Trans. Roy. Soc. (London) A238.229 (1939). 

2 *F. J. Corbato. Quarterly Progress Report, Solid State and Molecular Theory Group. M.LT., April 15. 1955. 
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. l 3 ? 'V hUn f ed ^ ,he ‘ e mU ‘ r ‘ 1 “ havc been imputed by Whirlwind. Many more will be done before 

the task Is done. This work Is being carried out by A. J. Freeman of the Solid State and Molecular Theory Group. 

262 N. EVALUATION OF TWO-CENTER MOLECULAR INTEGRALS 

V? to b * solved hy Whirlwind is the numerical evaluation of different types of two-center integral. 

T: y z£rjLX'° mic orbiu, “ ,or ,hr v “ iues of ,nternu< ,e “'- d — "• Th “ - -mg c“?£ 

for one .?"* *1“'"’ J* * v,lu,, “ 0 » ot one-electron Integral., has been completed. Routine, developed 

for one-elertron integral evaluation by F. J. Corb.tS of the M.I.T. Physics Department were used. ^ 

sut. ^ U ' C,UdBd *“ Ph - D - U « ta - * *• - - M.l.T. Solid 

263 C. FLIGHT PATH OF AN AIRCRAFT DURING PULL-UP 

xssz-jzr 

rdHr£S^ 

“* “ "**“ - 

computing lime using "Geor^ C* woukTbc^prohditUw^Therefore^^C^intc* 0 ^ 1 *! 0011 ' “ b * C ‘ Un ' “ PPttr ' n * ,hat lhe 
tic was written by Charles Block of the Instrumentation Laboratoi-y * ‘“’"l 

computation and the IBM Card Procram Culruint™- «. lh th id of * 8ma11 amoun t of hand 

Hon. of the CS program have be^flear.^^ ^ computational sec 

complete run of the flotations for the thirty sets of mtiaJ , onnfti remain * to b « d °ne before the 

ming errors in a cycling routine. ' 008 C ° n made *° remove some program* 

B.N. APPLICATION OF AUGMENTED PLANE WAVE METHOD TO CHROMIUM CRYSTAL 

Problem^;^ In Summary Report No. 37 under 

for any energy and angular momentum appear on the scope A routi "h" th *‘ lnteKra<ed " ave functions 
print, potentials and 22 for potentials written a. a ItaeTVnmM «' , " " rit,en " hich calculates and 

of the radius. The roulm P e ha. been applied to the potent^ or h e *P®"cntlaIs multiplied by powers 

of the Solid State and MoleeuUr Theory Gro' p’by «■ * Saffren 
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1. PUBLICATIONS 


Project Whirlwind technical reports and memorandums are routinely distributed to only a restricted 
group known to have a particular interest in the Project, and to ASTIA (Armed Services Technical Information 
Agency) Document Service Center. Knott Building. Dayton. Ohio. Requests for copies of individual reports should 
be made to ASTIA. 

The following is a list of memoranda published by the Scientific and Engineering Computations Group 
during the past quarter. 

No. 


DCL-41 

An ERA 1103-WWI Translation Program 

1-4 55 

J. M. Frankovich 

DCL-47 

Payroll Demonstration Routine 

1-13-55 

B. Rlskin 

DCL-48 

Automatic Scope Output Requests 

1-2555 

A. Siegel and S. Best 

DCL-49-1 

A Proposed Translation Program for the 
Numerically Controlled Milling Machine 

1-28-55 

A. Siegel 

DCL-57 

Calendar Demonstration Routine 

2-28-55 

R. J. Hamlin and E. Raiffu 

DCL-58 

Program for Solving Secular Equations 

3-15-55 

F. J. Corbato 

DCL-61 

Report on the WWI-ERA 1103 Input 

Translation Program (Published in the ERA 

1103 Newsletter) 

3-15-55 

J. M. Frankovich 

DC L - 6 2 

Report 18 January Meeting at MIT of WWI-1103 

3-15-55 

J. M. Frankovich 

DCL-63 

Annotated Prints of CS Flexo Program Tapes 

3-10-55 

J. M. Frankovich 


2. VISITORS 


Tours of the WWI Installation Include a showing of the film "Making Electrons Count.” a computer demon¬ 
stration. and an informal discussion of the major computer components. During the past quarter 13 groups visited 
the computer installation. Included In these groups were: 


January 8 
January 14 
January 21 
February 4 
February 14 
February 23 
March 22 


Professor Hansen's Class. Civil Engineering Department. M.I.T. 

Professor P. J. Rulon's Class. Mathematics Dept.. Harvard 

Reynolds Metal Company 

Atlantic Gelatin Company 

Union Carbide and Carbon Company 

American and New York Tel and Tel 

Air Force Reserve Research and Development Group 


The procedure of bolding Open House at the Digital Computer Laboratory on the first Tuesday of each 
month has continued during this period. Three groups totalling 65 persons visited the Laboratory at the Open 
House demonstrations. These persons represented members and friends of the M.I.T. community. the Retina 
Foundation, the University of Illinois and Harvard University. 

During the past quarter there were also 15 individuals who made tours of the computer Installation at 
different times. These Individuals represented French Morocco Realls Elect., Raytheon Mfg Co., Electronics 
Corp., A. O. Smith Corp., Dupont. Rensselaer Polytechnic Institute, University of California at Los Angeles, 
Harvard. Yale, and M.I.T. 
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